{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculate cross correlation functions from raw noise data in SAC\n",
    "\n",
    "In this notebook, we show some of the key steps in NoisePy for computing cross correlation functions in order for you better understand the unerlying processes of NoisePy script of S1. The input data are daily noise data recorded by two stations and here we show examples in SAC format. \n",
    "\n",
    "The steps to compuate cross correlation functions are:    \n",
    "* Load SAC data and the header info into memory\n",
    "* Break the continous data into small segments with overlaps\n",
    "* Perform Fourier Transform to convert signals into frequency-domain\n",
    "* Calculate cross correlation functions between the small time segments and choose to stack (substack) the cross correlation function of each segment and return to time domain\n",
    "* Save the function in an ASDF file\n",
    "\n",
    "More details on the descriptions of data processing, parameters for different cross correlation method and performance of NoisePy can be found in the online [documentations](https://noise-python.readthedocs.io/en/latest/) and our paper.\n",
    "\n",
    "`Jiang, C. and Denolle, M. NoisePy: a new high-performance python tool for seismic ambient noise seismology. In prep for _Seismological Research Letter_.`\n",
    "\n",
    "\n",
    "\n",
    "Chengxin Jiang & Marine Denolle\n",
    "\n",
    "Department of Earth and Planetary Science\n",
    "\n",
    "Harvard University\n",
    "\n",
    "November 2019"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Building env for NoisePy\n",
    "\n",
    "Before running this notebook, make sure that you have created and activated the conda env made for NoisePy. If not, you can create one using command lines below (note that jupyter is installed with the command lines here in order to run this notebook). \n",
    "\n",
    "```python\n",
    "conda create -n noisepy -c conda-forge python=3.7.3 numpy=1.16.2 numba pandas pycwt jupyter mpi4py=3.0.1\n",
    "conda activate noisepy\n",
    "pip install obspy pyasdf \n",
    "```\n",
    "\n",
    "Then you need to activate this notebook with the newly built NoisePy env by invoking the jupyter with the following command line.\n",
    "\n",
    "```python\n",
    "jupyter notebook\n",
    "```\n",
    "\n",
    "Now we can begin to load the modules needed for this practise. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import glob\n",
    "import obspy\n",
    "import scipy\n",
    "import pyasdf\n",
    "import numpy as np\n",
    "import noise_module\n",
    "import matplotlib.pyplot as plt "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Setup basic parameters\n",
    "\n",
    "The fundamental step is to setup the parameters used for cross correlation. As you can find from section below, there are many parameters needed for the computation, which are associated with the input data, some details controlling the targeted processing procedures and some tuning parameters. Some brief descriptions of the parameters are included here following the definination, but note that more details on this can be found in documentations.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "sfiles = glob.glob('./*.sac')      # find sac files\n",
    "if len(sfiles)<2:\n",
    "    raise ValueError('Abort! At least 2 sac files are needed!')\n",
    "outpath = './'                     # output dir\n",
    "\n",
    "# parameters of fft_cc\n",
    "cc_len    = 1800                    # window length (sec) to cut daily data into small segments\n",
    "step      = 450                     # overlapping (sec) between the sliding window\n",
    "smooth_N  = 10                      # number of points to be smoothed for running-mean average (time-domain)\n",
    "dt        = 0.05                    # sampling time intervals of the data: in real case it reads from data directly\n",
    "samp_freq = int(1/dt)               # sampling rate \n",
    "inc_hours = 24                      # basic length (hour) of the continous noise data        \n",
    "freqmin   = 0.1                     # frequency range\n",
    "freqmax   = 8            \n",
    "\n",
    "freq_norm   = 'rma'                 # rma-> running mean average for frequency-domain normalization\n",
    "time_norm   = 'no'                  # no-> no time-domain normalization; other options are 'rma' for running-mean and 'one-bit'\n",
    "cc_method   = 'xcorr'               # xcorr-> pure cross correlation; other option is 'decon'\n",
    "substack       = False              # sub-stack daily cross-correlation or not\n",
    "substack_len   = cc_len             # how long to stack over: need to be multiples of cc_len\n",
    "smoothspect_N  = 10                 # number of points to be smoothed for running-mean average (freq-domain)\n",
    "\n",
    "# cross-correlation parameters\n",
    "maxlag       = 100                  # time lag (sec) for the cross correlation functions\n",
    "max_over_std = 10                   # amplitude therahold to remove segments of spurious phases \n",
    "\n",
    "# group parameters into a dict\n",
    "fc_para={'samp_freq':samp_freq,'dt':dt,'cc_len':cc_len,'step':step,'freq_norm':freq_norm,'time_norm':time_norm,\\\n",
    "    'cc_method':cc_method,'maxlag':maxlag,'max_over_std':max_over_std,'inc_hours':inc_hours,'smooth_N':smooth_N,\\\n",
    "    'freqmin':freqmin,'freqmax':freqmax,'smoothspect_N':smoothspect_N,'substack':substack,\\\n",
    "    'substack_len':substack_len}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Load source data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read source and some meta info\n",
    "tr_source = obspy.read(sfiles[0])\n",
    "slon = tr_source[0].stats.sac['stlo']\n",
    "slat = tr_source[0].stats.sac['stla']\n",
    "\n",
    "# cut source traces into small segments and make statistics\n",
    "trace_stdS,dataS_t,dataS = noise_module.cut_trace_make_statis(fc_para,tr_source)\n",
    "# do fft to freq-domain\n",
    "source_white = noise_module.noise_processing(fc_para,dataS)\n",
    "source_white = np.conjugate(source_white)\n",
    "# num of frequency data\n",
    "Nfft = source_white.shape[1];Nfft2 = Nfft//2\n",
    "\n",
    "# find the right index of good signals\n",
    "sou_ind = np.where((trace_stdS<max_over_std)&(trace_stdS>0)&(np.isnan(trace_stdS)==0))[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Load receiver data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read receiver and some meta info\n",
    "tr_receiver = obspy.read(sfiles[1])\n",
    "rlon = tr_receiver[0].stats.sac['stlo']\n",
    "rlat = tr_receiver[0].stats.sac['stla']\n",
    "\n",
    "# work out distance between source and receiver\n",
    "dist,azi,baz = obspy.geodetics.base.gps2dist_azimuth(slat,slon,rlat,rlon)\n",
    "\n",
    "# cut source traces into small segments and make statistics\n",
    "trace_stdR,dataR_t,dataR = noise_module.cut_trace_make_statis(fc_para,tr_receiver)\n",
    "\n",
    "# do fft to freq-domain\n",
    "receiver_white = noise_module.noise_processing(fc_para,dataR)\n",
    "\n",
    "# find the right index of good signals\n",
    "rec_ind = np.where((trace_stdR<max_over_std)&(trace_stdR>0)&(np.isnan(trace_stdR)==0))[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Perform cross correlation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-0.00122514  0.00055406  0.00109561 ..., -0.00040456 -0.00018136\n",
      "  0.00081051] 4001\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEWCAYAAACT7WsrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd5wU5f3A8c/37jia9CZSPKSIiBXEjiiKoDEYo4klii0ajWn+EoMlqFiCxmhiNEajxhKNGktERVFRNGoop9KlgwJSjiJF6nHf3x/zzN7s7syW270CfN+v171ud+aZmWdnZ+c7T5lnRFUxxhhjwhTUdgaMMcbUXRYkjDHGRLIgYYwxJpIFCWOMMZEsSBhjjIlkQcIYY0wkCxIGEblYRD7KYfk3RWRYPvPk1ttQRF4TkfUi8u98rz/NtmeKyIBqWO94Ebk83+s1NWtP+h4tSJisiMgtIvLP4DRVHaKqT1bD5s4G2gGtVPWcalg/ACLyhIjcHpymqgeq6vjq2mZV7EonJpfXdSJSPzDtchH5ImFaKxFZJSKDRWSAiCyNWNfl7vUAEVEReSUhzSFu+vhq/FgZcb+RHSKyKfD3TWC+ish0ESkITLtdRJ5wr/+WsOwmEdnslusfWMcqESkKrKOem5bXm98sSGQg+EXUZeIpSDdtF7IvMFdVy2s7IyZzIlICHA8o8F1/uqo+CiwDRgSS/wkYo6pvZbGJMuBoEWkVmDYMmFvFLFeH51V1r8Bf84T5+wDnhi2oqj9JWHYv4GXgfeDjQNJ1wJDA+yFuWl7tqiePvBCRTiLysoiUicgaEXnATb9YRD4WkftEZA1wi4gUiMhNIvKli9ZPiUgzl76BiPzTreMbEZksIu0C61ooIhtFZJGIXBCRl0IRuUFEFri0n4pIJzfvGLfO9e7/MYHlxovIHSLyMbAZ2C9iWjMReUxElovIMnflUhiRlz+LyBIR2eDycbybPhi4Afihu7qZGsiDf6WXaj+VuCugYSLylYisFpEbI/JwK97JxN/WZZJQigmsryiQj9vcd7dRRN4WkdaB9MeJyCfuO1rivpsrgAuA69x2XnNpF4vIye51fRH5k4h87f7+JO5qWNzVr4j8n/u8y0XkkohDztdVRCa5/fuqiLQM5PGoQB6niqvyEpE78E68D7h8PiAit4rIX9z8eiLyrYj8wb1vKCJb/XVHrdfNizw23D76SETuEa9ksEhEgiemMBcBE4An8E7eQZcDV4vIoSJyKjAQ+FWa9SXaDvwHd5J1ef0h8EyqhUTk3yKywv2OPhSRAwPznhCRB0XkDXfsTBSRroH5p4jIbLfsA4BkmedEdwO3SgYXoCJyFXAicJ6q7gzMehpvX/suAp7KMV/JVHWP/AMKganAfUBjoAFwnJt3MVAO/AwoAhoClwLzgf0AP7I/7dJfCbwGNHLr7QM0devdAOzv0rUHDozIz2+A6cD+eAfgIUAroCXe1cGFLi/nufet3HLjga+AA938ehHTXgEednlqC0wCrgx83o8CefmR23YR8H/ACqCBm3cL8M+EvI8HLnevU+2nEryry7+7fXoIsA04IGKfxG0r5L2/vqJAPhYAPdz6xwOj3Lx9gY1u/9Vzn+9QN+8J4PaEbS8GTnavR+Kd9NoCbYBPgNvcvAF4x8pIt97T8AJzi4jPNB7varq3+y5e8j8T0AFY49ZRAJzi3rdJ3M/u/UnAdPf6GPfZJwbmTc1wvemOjR3Aj/GO7auArwFJ8duaD1yN9zvYAbRLmP8z4DNgEXBmYPoAYGnEPrs8mMZ9Xv+zngaMxQtA41Pk61KgCVAfrwQzJTDvCbdP+uEd988Az7l5rfGOnbPdd/wr951fHrGdW0j4jSTMV6A78Gngc90OPBGS9ghgE+7clLCO3sBKoDnQwr3uDWhez5X5XNmu9AccjVdsLQqZdzHwVcK0ccDVgff7ux9AkTv4PgEOTlimMfAN8H2gYZr8zAGGhky/EJiUMO1/wMXu9XhgZML8uGl49frbgnnAO1m+H/i8H6XI2zrgEPc66QeQ8CNOtZ9K3MHdMTB/EnBuxHbjthXy3l9fMEjcFJh/NfCWe3098ErEdp4gdZBYAJwWmHcqsNi9HgBsCR5HwCrgqIhtjccFLve+F96VcSHwW1xADcwfCwxL3M/ufUNgK17AG45XyluKF5xvBe536SLXm+GxMT8wr5Hb53tHfL7j3Pfd2r2fDfwqIY0AExO/D7cvK/B+M8G/2AmZQCAB5uEdX8/hlQZTBomEbTV3n6NZ4Bh4NDD/NGC2e30RMCEh/0tJHSS2J3yG9wPzFejmtvElUExIkMC7QFwMXBuyDX8dj+JdpP4E7+KrG3kOEntydVMn4EuNru9ekvB+H7wv1Pcl3omvHV6xbyzwnKuOuFtE6qnqt3jF4J8Ay11RtmeK/CwImZ64XX/bHVLkNXHavnhXQMtddcM3eFeObcMyIiK/Fq+Bcb1L2wzvaioTqfaTb0Xg9Wa8k1q+RK07av9mIuwz7RN4vybhOEr3mYLfzZd4301rvO/pHP87cvv+OLwSaBJV3QKUAicA/YEP8C5WjnXTPnBJU603k2Mjtk9VdbN7GfX5hgFvq+pq9/5ZEqqc1DvLfQHMDFn+a1VtHvwDonrePQ1cg1cV80pEGiBWnTtKvOrcDXgnX4g/rqOOnX0IfGcu/2G/uaAXEj7HiYkJVHUMXrC5MiS/AvwTKFXVe1Ns5ym8IFY9VU14P9491RKgs4gURQQKTXj/Nd4PytcZ7wpnpVv+Vrw6xhJgDF7J4DFVHQuMFZGGeFcLf8erWw7LT1dgRprt+tsONvQl5jVx2hK8q8XWKYIiAOK1P1yHV1c8U1UrRGQdlXWwYdtKld/YfgI6plk2nW/xrmR9e2ex7BK8qoQwmX4m/6TW2U2rqk6B153xrrxXuzw+rao/ziKfH+BVLR0GTHbvT8X7rB+6NJHrFZH2ZHhspOOO8R8AhSLin3DrA81F5BBVnZrL+kM8jVe19ZSqbvbOq5HOB4YCJ+MFiGZ4JeRM2haWE/jO3Am8U3TyrNwI/Mv9Bd2EVyrom2b5/+IFe8ULpl1TJ8/enlySmIT35Y8SkcbiNT4fmyL9v4BfiUgXEdkLuBOvB0O5iJwoIge5BrQNeD/6ChFpJyJDRaQx3g9xE15xOsyjwG0i0l08B4vXe2MM0ENEzheRIhH5IV4VxeuZflBVXQ68DfxRRJqK17jcVUROCEneBO+kXgYUicgIvPYV30qgRKJ7TEXup0zzm8IUoL+IdBavMfz6LJZ9BjhZRH7g9mMrETnUzVuJ14YS5V/ATSLSRryG8BF4V3lV9SMR6SUijfDaMl5Ur0Hyn8AZInKqu/JtIF7DuB9cw/L5Ad5V5CxV3Y6rkgIWqWqZSxO53iyPjXTOBHbiHZ+Hur8D8E5kF6VYrkpUdRFeiSm080OCJni/wTV4Fxp3ZrGpN4ADReQs19D8c7K7QImkXjfrGQRKW+J1mLgO+L6qbkizvAJnAN91r/Nujw0S7kd5Bl60/gqv2PfDFIs8jnfl8iFeg9tWvAY48A6YF/ECxBd4P9yn8fbvtXhXnWvxDuirItZ/L/AC3g92A/AYXj3xGuA7eA3Ia/AOnu8EivOZugiv7nMW3hXUi4RXY4zFK6XMxasK2Up80dq/qW2NiHwWsnyq/ZQTVX0HeB6Yhtfol02g/AqvDvj/8L6LKXgN5+Dt616uuuU/IYvfjletMw2vc8FnblpVPY1XB74Cr8PEz10el+Bd7d6AF6SX4HVo8H+nfwbOdr2M7nfTPsFrm/BLDbPw9rn/PpP1ZnpspDMM+IeqfqWqK/w/4AHggkx68mRLVT9S1UxKdU/hHc/L8D7nhCy2sRo4BxiF9xvsTnxX1DB+r7zgX2j1Ll6poWXg/Q143+n/QtaR1DtSVWeqaljVXV5INQUfY4wxu4E9tiRhjDEmPQsSxhhjIlmQMMYYE8mChDHGmEh56W0g3pg+f8a7a/RRVR2VML8+Xu+CPni9A36oqotF5BS8HgPFeHco/kZV33PL9MHrAdIQrxvoL9J18WrdurWWlJTk4yMZY8we49NPP12tqm3C5uUcJNy9AQ/ijQezFJgsIqNVdVYg2WXAOlXtJiLnAnfhdTddDZyhql+LSG+87pf+ncQP4Y0XMxEvSAwG3kyVl5KSEkpLS3P9SMYYs0cRkcRRHWLyUd3UD29sl4XuZp7n8PpkBw0F/OcNvAgMFBFR1c8DfZxnAg3FG3GzPdBUVSe40sNTeDfqGGOMqUH5CBIdiL/Zainx4wrFpXF33q7HG5Qs6PvAZ6q6zaUPPnwkbJ0AiMgVIlIqIqVlZWVhSYwxxlRRnWi4Fm9c97sIGegqHVV9RFX7qmrfNm1Cq9SMMcZUUT6CxDLiB7vq6KaFpnG35jfDa8DGjUvzCnCRqi4IpA8OBhe2TmOMMdUsH0FiMtDdDehWjPe0qNEJaUZTOYDV2cB7qqoi0hxv8KzhqhobC8UNOrZBvKdpCd7YMq/mIa/GGGOykHOQcG0M1+D1TPoCbxz1mSIyUkT859s+BrQSkfl4A94Nd9OvwRtgb4SITHF//iBYV+ONjDof7zkAKXs2GWOMyb/daoC/vn37qnWBNcaY7IjIp6oa+uyKOtFwbUxdsu7b7bwxbXltZ8OYOsGChDEJLn+qlJ8++xlrNm2r7awYU+ssSBiTYO7KjQAUpH4cpjF7BAsSxiSoqNh92umMyZUFCWMiWKgwxoKEMUnEVTNV7EY9/4ypKgsSxiTwu4VbjDDGgoQxkdQqnIyxIGFMIr+6yUoSxliQMCaSBQljLEgYk8S/O8Iaro2xIGFMEk34b8yezIKEMRF2p8EvjakqCxLGJPCrmyxGGGNBwphIFiSMsSBhTDJXlLCGa2MsSBiTTOP+GbNHsyBhTARruDYmT0FCRAaLyBwRmS8iw0Pm1xeR5938iSJS4qa3EpH3RWSTiDyQsMx4t87EZ18bU71i1U21mw1j6oKiXFcgIoXAg8ApwFJgsoiMVtVZgWSXAetUtZuInAvcBfwQ2Ar8Dujt/hJdoKr20GpTSyxKGJOPkkQ/YL6qLlTV7cBzwNCENEOBJ93rF4GBIiKq+q2qfoQXLIypEyrvuK7VbBhTJ+QjSHQAlgTeL3XTQtOoajmwHmiVwbr/4aqaficS/ixJEblCREpFpLSsrCz73BuTIHbHtQUJY+p0w/UFqnoQcLz7uzAskao+oqp9VbVvmzZtajSDZvdmQ4Ubk58gsQzoFHjf0U0LTSMiRUAzYE2qlarqMvd/I/AsXrWWMdUuVt1UUavZMKZOyEeQmAx0F5EuIlIMnAuMTkgzGhjmXp8NvKcp+heKSJGItHav6wHfAWbkIa/GZMxKEsbkoXeTqpaLyDXAWKAQeFxVZ4rISKBUVUcDjwFPi8h8YC1eIAFARBYDTYFiETkTGAR8CYx1AaIQeBf4e655NSYb1iZhTB6CBICqjgHGJEwbEXi9FTgnYtmSiNX2yUfejKkqCxLG1O2Ga2NqlVU3GWNBwpgkfm9ru0/CGAsSxkSysZuMsSBhTCQLEcZYkDAmkpUkjLEgYUwkixHGWJAwJonYUOHGxFiQMCaCVTcZY0HCmEgWIoyxIGFMpAorSRhjQcKYSBYjjLEgYUwUixHGWJAwJknl40stTBhjQcKYCBYjjLEgYUwkK0kYY0HCmEgWIoyxIGFMNIsSxliQMCZR5fMkLEoYY0HCmAgWI4zJU5AQkcEiMkdE5ovI8JD59UXkeTd/ooiUuOmtROR9EdkkIg8kLNNHRKa7Ze4X//LOmBpiJQlj8hAkRKQQeBAYAvQCzhORXgnJLgPWqWo34D7gLjd9K/A74Nchq34I+DHQ3f0NzjWvxmTDQoQx+SlJ9APmq+pCVd0OPAcMTUgzFHjSvX4RGCgioqrfqupHeMEiRkTaA01VdYJ6Q3E+BZyZh7wakzErSBiTnyDRAVgSeL/UTQtNo6rlwHqgVZp1Lk2zTgBE5AoRKRWR0rKysiyzbkwyv17Thgo3ZjdouFbVR1S1r6r2bdOmTW1nx+wGNOG/MXuyfASJZUCnwPuOblpoGhEpApoBa9Kss2OadRpTrazh2pj8BInJQHcR6SIixcC5wOiENKOBYe712cB7mqIsr6rLgQ0icpTr1XQR8Goe8mpMWpXVTbWaDWPqhKJcV6Cq5SJyDTAWKAQeV9WZIjISKFXV0cBjwNMiMh9YixdIABCRxUBToFhEzgQGqeos4GrgCaAh8Kb7M6bGWIwwJg9BAkBVxwBjEqaNCLzeCpwTsWxJxPRSoHc+8mdMNvw7cqzh2pjdoOHamHzzY4PFCGMsSBgTSa3CyRgLEsYk8qubKipqNx/G1AUWJIyJYOUIYyxIGBPChgo3xmdBwpgkGvfPmD2ZBQljIljDtTEWJIwJ4Vc31XI2jKkDLEgYE8GaJIyxIGFMklgXWIsSxuRnWA5jdgdjZ67gsY8WVd5xXbvZMaZOsCBhjHPl058C0KJRPW+ClSSMseomYxL5DdbWcG2MBQljklS46GCjwBpjQcKYJH6DtYUIYyxIGJNkpwsSVt1kjAUJY5L4o79adZMxFiSMSeKXJCxGGJOnICEig0VkjojMF5HhIfPri8jzbv5EESkJzLveTZ8jIqcGpi8WkekiMkVESvORT2MysdNvuE5olVj2zRZUldtfn8X8VRtrI2vG1Licg4SIFAIPAkOAXsB5ItIrIdllwDpV7QbcB9zllu0FnAscCAwG/urW5ztRVQ9V1b655tOYbAVLEjOWrefYUe9x11tzePSjRVz8j8m1lzFjalA+ShL9gPmqulBVtwPPAUMT0gwFnnSvXwQGioi46c+p6jZVXQTMd+szptYFG66//mYLAOPnrAKsKsrsOfIRJDoASwLvl7ppoWlUtRxYD7RKs6wCb4vIpyJyRdTGReQKESkVkdKysrKcPogxQd9s3s4to2eyvbwCcQM61XRw2Ly9vGY3aEyCutxwfZyqHo5XjfVTEekflkhVH1HVvqrat02bNjWbQ7Nbe/jDhTzxyWJenbLMDR5es8+YmLFsPb1GjOXN6ctrbJvGJMpHkFgGdAq87+imhaYRkSKgGbAm1bKq6v9fBbyCVUOZDCxZu5mpS77J6zorVGMjw9ZkSWLGsvUAvO+quIypDfkIEpOB7iLSRUSK8RqiRyekGQ0Mc6/PBt5TrxP6aOBc1/upC9AdmCQijUWkCYCINAYGATPykFezmzv+7vcZ+uDHeV9vLEjkfc3Rigq9n2d5yF19z0z8kjemWQnDVL+cR4FV1XIRuQYYCxQCj6vqTBEZCZSq6mjgMeBpEZkPrMULJLh0LwCzgHLgp6q6U0TaAa+4euAi4FlVfSvXvBpTFarE2iTmr9qUMm1FhbJxaznN/JFkc+BiRKxLbtCNr3jXTKcffHrO2zEmlbwMFa6qY4AxCdNGBF5vBc6JWPYO4I6EaQuBQ/KRN2PyQdInAeBP787l/vfmM/u2wTSoV5h+gRQKXGAKCxJB6zfvoEFxAfWLctueMWHqcsO1MXWCUlmS8ElE1Hjso0UAfLst915JfpBIV8V1yMi3Oe+RCTlvL99UlZv+M51pS/PbRmRqlgUJs8e4+plPeXVKYp+K9FSTSxJL120JTVtY4KXcWu4NAPVC6RJKhr/BJhc0tpdXMHrq12nHhVq1cSvfbNnhtp++JeSzr+reiXjDlnL+OeErfvToxNrOismBBQlT6yoqlEWrv62WdX+55ttYYBgzfQW/eG5KldYTVXJIVOyqfLbu2AnAwx8sAGC5uxnv/nHz+Pm/PuedWStTrqffHeP43X92j74adt/hrs2ChKl1fx0/nxPvGc+cFZmPh7Ro9bes2bQtbbohf/5vlQODT1FmL88sb8WFriThgoQkVBl9vd4LFuu37OCXz31O39vfTb/9XfUsW3lzidmFWZAwta70y3UALPtmc8bLnHjPeI6/+/206TZv35lyfvnOitgJPcq0Jeu5Y8wXGeXLDwrlO70zo3+e/PyrdSz7prKKasdO5T9Tvmb1pm3c9dZsVm7YGrnOYJBYvWkbI17dNUoYmZa+TN1mQcLUutgFZ5ZXnJu37+TW12bmtO3z/j6Bnr9L3bv6+dIlKeeH8Ycb90+Uv31pOseOei92VX3z6MoT/UPjF/Drf0+NXFfwLu/bX5/FU//7Muv81KZsCxITFq5hSp5viNzdfDC3LDaOWHWzIFEL1m/ZwewVG2o7G6G2l1ek7XKZb7mMi/SPjxfntO3Ji9fltHyUiQvXAiAJTd7+R9yxM/7Dbi+vYPLitbGBBKPsik/Ly/bhTec+MoEzq+GGyCjzV21kVYqSXE3asn0nC8pS34sDMOzxSTU2ErEFiQxMWrSWo38/jpLhb1C+syLn9Z37yAQG/+m/WS0z8+v1DH9pGhVVOEts2b6Tdd9uD523s0LZXl75mXrc9Cbn/z19d8qyjdtYui7z6qFUNm51vXgi5m/dsZNuN4yJ7JkUdhJaG/F5o8xdmbrN4dJjuyRN+/yrdXEnl2+3lbPe9Ui6663ZQHKVyyufR/euOudv/+PEe8YnTQ9+vIJdqAqnKkF/YQYnSIBP5q/mofEL3HaUa5+fwoSFa7LfIHDyvR9y1O/HVWnZXCxZuznpmP7ps58x8I8fZHye2bpjZ5XOCdmwIJGBHzz8P5av904G28pzDxJfLM++FHHpE5N5bvISVm1M31ib6Ht//ZjDbnsndN5Fj0+kx01vsrNC2eBO1hMXreXfpUtiJ+9Ey9dv4Yg73uW4u9K3CWTCv5qPuuIs27iN8grl7rfmhM7fsLWce9+ew3UvVlbZbE/4noLrXrT6W065N/6HOOi+D1PmsW3T+knTvvfXTzj53g9i70/4w/hYV9dsbXHtImHHV0VckEgfJVSVr9bEB/ANW3dQMvwNXpv6dZXyl4lPFqymZPgbzPMDrsb9y8hJf/wgfSLg/EcnxgLxjp3Ky58vy6mrbW2U0M588OOkThUfzPVGst6ZYYTt+bu3uKma26gsSCRQVe56a3ZkT5tsjqWtO3ayIw8lD4CVG5KDwz1j5zAmgxFCZyd8ll889zmXP+kVVT+e7119jXh1Bgff8nYszW9enMZNEV0w565Mf7X36pRl3Pb6rLTpgqr6O73vHe8u5xdKl8amlVdUxLVXBE8CJ94znnmrNnH9y9Mz3kZFxI92w9bKoLA6pLdV4k14UaYtXR97/daMxO+0ctsFCUWJsMD68IcL6f+H9ykZ/ga3vT6L7eUVfLl6s5u3gB07K/i/F6by3KSvWLJ2M3NWbMz4Cj6V191YUhMWrXW5zu4xsLO+rloVrP/dBANovzve5YmPF8Xel++sqFLJd/ycVSHfR7RXpyzjwseSg9XGrTt4+TPv+FRVKiqUNa60G/wOY09FTNhnM5atj/s8Qc9O/Crj/FVFXobl2J1s2FLOQ+MX8Nykr/h8xKCkK9J0vli+gXWbt3NM19b0/N1b9Ny7CW/9MnSUc8ALJEUFEhvMLZ2jfj+O1392HL07NOOB9+cDsHhUZuP3rNqwlbZNG/DqlOSryf+EVIOURZRaUp32Hv9oEV9/s4VH3Z3HN51+AI99tIh+XVpSXFRAz72bRi575dOfcl6/zvz+rINSf5AEYb2TKiri2yumL1uflObLNZmfNFIV6TdtK6d08dqk6SXD38h4/UF/ende5LzE6qaw9iP/rm//dY92e9GrfbPYtMmL1/LSZ0t56bOlcctlehwlUlUmL14XO9n5ecymuklVOe3+7KpgAa54qpSfndQd8Kr2lq/fwjGj3kMVbnltFhe7asK7x87hkQ8XMumGgbRt2iDtej+cW8ZFj0+Kvc903/glg0/mr+aYbq1j0294ZQavTf2a7m2bMGnx2rgLqAqFwoTvNXHffecvHwFwcq92VOTnujNjFiScVRu20qRB5aBs/o8v8erjm83bKSqQ2Lg8677dzrJvttC7g/cjHPJn70D3D6rgVfyqjVt59L+VP2BVpefv3uLUA9vx8IXJT2hdv3kHxUUFNCyOH5Nn9NSvY9vLRr87x/HutSeEzsukGsMXlbRs4zZGJpQeZn69gdvfqOw+6u+Xb7eV07h+EdvK40/w/5r0FacdtDfHd2/Dq1OW0apxffZt1SjldsOu1tdtjm+T+PTL5AbqTIv0AAtT3Oz3q+enpL05LhuJ2fLff/rlWjZuja/OCotdicE9WJidu3IT9YuqXoFQMvwNLj+uCzd9pxc3vDI97iq2uRvU0G+s97O2ZcdO3py+nC9WbGTM9OW8/rPjYr+fTxaspmmDevRo1yRpW3e9NZsr++9H80bFjJ+zivvencfLVx0Tu6sd4O1ZK3nb7XsRGD+nLG7/vTB5Cef07cgjHy4EYPWm7RkFiWCAAK/9oFPLRmmX853/6EQaFRcya+RggFjb1RkPfMT+CZ91Z4XGfSaIfm5Jvqp4s2HVTU6/O8dxwIi3GDfbO+D8ryjxyzvurvfjel6c9dAnsSj/3uzKE0Ww4fTbbeW8PXMF/e4YFztYofIHPnZm+AnmkJFvc8CIt2KBxxd2rqyoUCZm0HA3LOHg920MqUufv2pTaB37ivXJPUFUNbTksXhN/Mn16mc+ZeLCNRx481g+mFvGXW8mtzNc+NgkyjZu4xfPTeFHj02MnfCXrttCyfA3WL85vq0kLHgkDhceVvWVTRvwy59FNzhn0hslG3MSGtEV7/v9/kP/480ZK+Lm/Xnc3Njr9Vt2hFY/BQ/h7eUV/PL5zG8u3LJ9J+/NXomqxqrT/FJiYjXHN4HvZezMFRweaAe76pnPuH/cPOav2sS4L1bR48Y3mbLkG87/+0S+85ePQu8TeWj8Ag4d6a3jty9NY+qSb1LeTyJIUoC97qVpjJleuc92Vmhlm0kWgvfkqGpGN35u3r6TRau/5Z8TvmTiosqSZuL3G9ajrS71YrOSRIJrX/AaPzduLeerNZuTggRUlg5WrN8aN5zEpU+Uxl4HfyAH3jw2dFtR9dyJEhu6P5hbltTA+ehHC7lzzGx+0LcjN57Wi807ymnVuD7FCVeNy9J0sQxatXEbl/xjEqs2bqN5w3q8es1xjJm+nN+8OC2WZvP2cgpEeKF0CSNeTb5n4ZpnP497P2b6itsD7c8AACAASURBVNiP9pP5q/lqbfgV+umBqoe/uaEtfK9Pj68u27Q1u8bi/j3a8OHcMn54RKfYjXx1mapGlnoefL9y3xxy69vc/f2Dk9K8OuVrDmhfWc23ZG3qY2DHzgremLacoYfuwwEjwu8hGfdFdMlJBP6cosrso/llbN9ZEXexle7GSL/kkar6d8uOnTzoqmCD/B5n4F3JA7z1y+Mjqz6jehYtX7+F9s0acu87c/nLe/N54cqj6delJcNfmsZzk5eEVkmF9VZLNOCe8SwedTpL1lbWWvjnhk3byrM+vvPNgkQK/f/wPj87qVvovNkrNmTdjTVRsD55646dGQ8tPXvFxqTGaL/R74XSpXy1djMTFq7lrMM6VKlaKsjvefQlMG/lRq5+5rO4+b1GeAHw+O6tExdN6+FAqSpRsBdX8EoQvC7JQaOz7LHzlSvdJAbQqlpYVj3jTvm+3baToQ9kdt/Aa9OS98X/Fq7hgix6/jz8wQLueXturMdVmMueLI2cl65DQLb34Rx557uxjhuvTf2anw3sHpk27CIosUoTYPk3W1m2bkvo57jn7blJ0wCO/v17TBlxCn95zwtEj/53Ib9/8ws+d4MrBtuCsrW9vCKuwfvyJ0u55+xD6P+HzKqXHvlwAVf071rl7adi1U3A+7Oj71z0D4hEiQHi9ZAfZzrBq+Wev3uLt2euYNO2clZv2sbZD32S8XomL17LjsAPb/4q76T18ufLktoIgrK9W/mUFN1E/ztvdVbrykWuPcYWuwbrXMd0qimTFq9lVobdpldvCr8/JJOuuX5VlX+SzKb3VzaCvdAyEezZ98d3wk/gqdz6WvhvIDFAVFQo4+es4tMvkzsh+O4fV3k+eHvWyliAgPAqzUz1uOnN2HEJ3oXQgtWZV2PeOWY2y9dnXkuQDStJQKyXUC4Sq1UysSDhCvSKpz+t0rbvfXsu/wu0R4R1xQyT693KtSWxZGEqVeUeHF/XG8Zw69DeecxN9fBvostF2IXGfjeMCUkZ7/GIbqjVYefO7EpcT3yymOuHHJD3fEi2t8zXZX379tXS0uhicJSzH/pkl6ibNsbsOeoVStLwLakcUdKCf//kmCptS0Q+VdXkLpZYdZMxxtRJ2QQIyK4be1brzcdKRGSwiMwRkfkiMjxkfn0Red7NnygiJYF517vpc0Tk1EzXmU+7T1nKGLOnKkq8Iy9Pcg4SIlIIPAgMAXoB54lIr4RklwHrVLUbcB9wl1u2F3AucCAwGPiriBRmuM68Ka9LnZKNMaYK6nJJoh8wX1UXqup24DlgaEKaocCT7vWLwEDxbpMdCjynqttUdREw360vk3XmTT5GdjXGmNoUdk9XPuQjSHQAgk9lWeqmhaZR1XJgPdAqxbKZrBMAEblCREpFpLSsrKxKHyBfg/AZY0xtKazDJYlapaqPqGpfVe3bpk2bKq2jPMsGImOMqWsa16+eOxryESSWAZ0C7zu6aaFpRKQIaAasSbFsJuvMG2uTMLuS35y6f21nwdRBPz0xfHSIXOUjSEwGuotIFxEpxmuIHp2QZjQwzL0+G3hPvRs0RgPnut5PXYDuwKQM15k3RbvS477MHq9ji4ZceNS+tZ0NU8c0bVhHSxKujeEaYCzwBfCCqs4UkZEi8l2X7DGglYjMB64FhrtlZwIvALOAt4CfqurOqHXmmtcoPziiU/pENejgjrmNt2Sq37OXH1lr6y8skGprpDS7rrrccI2qjlHVHqraVVXvcNNGqOpo93qrqp6jqt1UtZ+qLgwse4dbbn9VfTPVOqtLdTX4ZOucPh0BGNx775zX1WffFjmvY3fxyIV98r7ObJ5FURXBB9YkKhSJPWMjTLrnRZx6YLukaZ2zeFZCTTjrsNB+KruVqSMGVXnZrm0aJ00rKqieJuZdvuE6H3L5wZ92UPwJvWOLhlVeV8vGxUDlQ1uq6soT9uOlq6p2e36m6koVXZsm8c+ennHrqUlp/P2aT7XZjFVQIAw7uiRy/rEpAgyEH1+XHFtS5SfTVYd8jdBbnZ7JsTTZrFG99ImAwzs3T5oWVmqo0yWJXd0Zh+zDzWck36vXNuEEBPDaNcfFXrdrWp8D94mvGvrzuYeFbqN3B2/s+m5t9wKSgwtUDqGc+F2ffEDbFLn3/HZwz9jrE/dPnz4T/XuE9xa79weHhJ6M06lXKEy/pepXT2EeubAPjQJP7gv7mVTHj2e/1o1zesJbLgpFKCiQyNJiYYHwq5N7MHxIz9D5325PHhH29IPbA/Dkpf04pFPySQng2lN6JE1b9PvTUua1R7u9Us6Pkumw+cUJj/1t0ageX4wczFmHdYg9Kc+X+EQ430MXHF6lPAaD8aXuMan59twVRzH00ORSVeeW8SWJ4qICGhdnts+yZUEC6NC8IRccuS/HdG3Fj47qDMBJPdvy0W9Piks3pPfeHBRoLygQSXoSWIfm8SWJn5/UjZvP6MXonx7H4lGn89cLDuf47q05vHP8D/z5K46K9bKql3Dgn3V4x9jrf/34qNDPEDwPHrVfq1QfN1TffVtQXFQQdzX51KX9QtOedXjHpB/nwjtP45iurfhh3+j2nV8M7B73iNigdk2TA/KwoysbZ/9x8RGhyzVtWC+uujDsrtPCAuGswztwXLfWse83VwUFwi9PTj5pplPVk2aQH/SiQl9RgfCLk7tzybElofMTh3VfPOp02jbxHul5Qo82PH9F+DH285DnOIQ9Ojbo8cD3lkmvrE9vOpnbzuzN/nuHn9CDpow4hUk3DuSD3wyI5W3YMSU0LC7k3h8eypQRg7jhtMpAWVxUwG1DD0xaz3FZPAvlk+EnhU4fEXKReU6fjlw1oCvHpSjZhQXkQb0qqwOP2q8V5x/ZmasGeM+K2Kt+EQ9dcDijvh//HPi5tw+hqNCqm6pVcVEBz/74qNgJ+aSebWMnzQfPj7/SeOIS78C/9NguFEQ8m7aoQHjhyqO5dtD+XBJI16NdE56+7Mi4QPDzgd05cr9WdHcnkP3aNOYX7qC/6/sH0dBdVTUuLuTorq34/HenJOU/7Lf63+tOjBXbbzq9cgjhV396bOx1p5ZeUHt0WF/m3j4kcv8knsQLCoTDAsXgggLh2R8fxV1nH8zUEYN49sfJRXH/hBL29LSJN5zMX86LL4WVtG7MS1cdw8ihB3Jiz7a8fHVyFZoAhYExa8L2Q4EI9/7gUP55+ZFs3ZGfGycLJPo5xKkU5qHe2D+WooZh8INIvSpuK1XJa9INA2Ov/SrHl646OjbthSuPjku/V6Dv/k9P7MYRJanbylrtVZ8Lj9qX77iSTZTmjerRvFExzRsVs2+r5Pp5X2L39lN6xZfgz+vXiSYN6oXWJJx8QLvYb93XLuH52Cfu34Z67vjr1LIhffZtwXcObk+TBkX8dkhPfju4J/9MUS2VuKsLBB65KH4w1nqFBVzjurdWqDLkoPZx+/X6iBJjvtjzJBIc3rkFH/7mxNjJM8yA/dvGrrg3by/n62+28M8J3vN+/ZP/4fu2oF+XlpHrCD71yz9Ozu/XmYM7NOegjs04oUcbrhrQlQb1Ctm8vZx+JS251V0FtWhczDUndot7DoYgfOfg9pzUs7KqqVPLRnzv0A48X7qEc/p04rDOLVi6bnPc1ctDF/Rh9aZtNG8UXW/foXlDXr3mWPre/m7c9FeuPpaS4W8kpW/WqB7HdG1N97Z7MW9V8oNT+kacKM44ZB8GHtA29rS7c/p2Yq/6RbFqlcTSF3iDMwbbR9L1QTivXyde/DS7h96E8UqRydN77t2Ee845JPbc80QVeWjM8EtOUZ/1e67Rt6BAuPG0A3h16jJmLKt8zkTXNo2TnmUStn7wSpjBYfTbBk6Sb/z8eAD67NuSD39zIh1aNIws3TRpkHyqeffa/nw4dzVn9+3Iwbe8nZC+Ho9c2CfyGSuJ+35gz7bcP25eUlXrqQfuzd1vVT5Hfe9mDfhk+Em8M2slR3dtRXdX/XvJsV2YvnQ9L3++jH//5Gi6tdmLFo2L2bK98nd64VH7xk7qfon0H5dUlrb/e114KSPogPZN45754Qf6F648mh88/L+4auOSQOcEP3D7VavBQH5F//3SbjcXVpII0blVo9BidNiPslFxEbefWVn0a71XfR6+sA8P/yh1j5pzA91u/WoBEYlVZ4lIrF62UXERL/zk6LjnFCdexYrAA+cfHlc1BXDbmb0Z/+sBNGtUjz77tkiq32zTpD4DUrRhzL9jCP+97kRa75VcHZTOO9eewOJRp/OTE+Ifqxg8Tx7TtRW/HlRZbdOouCiujSFK8Af0+7MqSybVNchZIpHwE/7PB3and4dmSY3AJ+7vte/4nSQO3KcpZx1etR48hSlKEmf36cjAAyqrK37cfz96tPWqbkYOPZBptwyKe3zo+78ekLQOv6RyXr/O/Cui6gmIqxLq3KoRhQVeW8ktIVflfk79q/7BB+5Nt7ZNuPS4LjSNqIIcdODeST2x/BJ24uNyD+nUnMWjTk+qvunaZi9GX+OVnP3fzD7NGzLsmBJ6tGsS9zu/86yDGPvL/hxR0pIWrrNDw+LCWFvGVQO6IiLMvX0II7+b3cOZ7vheb7q13Yt7zokvRfvnelVl8ajTudL9VibdMJDXXRAGr43mhtN68twVXkktGMjTVfnlykoSGfBvUvHrbdM59cD0XVibNyrOqTfJkN7tefD9BRzTtRWfLFgTma64qICS1tHF8cTiszetfuyRkfms56w8lr0fa8+9m/BsSBtLr/ZNKf1yXcquyf4JUhVOCdThBk+cPzqqM/+c8FXKRtC7v38w1700LYtPUbmdHxzRKelxmlG9vo7p2pr355TFTtDtmjagecPw0tsbP/c6R0y4fiBrvt3G6ffHl0r86o2w2qRRZx2UNM0PCY2Li2jaoF4sUH34mxPpHNGVds7tg6lXUEBBgfDJ8JPiqkev7L8f736xMnQ5gIuP7cLAA9qxfWdF7IrfDzz+eo7vEX+S/8clR1ASUm308IV9Wb95B4eM9EoavzqlB2cd3iH0uI2ydzMv7am9Uv8uG9QrDG0LefnqY/j8q2/Yx7U3VqXn1QVH7ssFR+7Lhq074qb7J/jE6422IZ8v+AzrxGru6mRBIgPHdWvNfT88hCG9U9eT1iT/avW1qV/zyYI19NqnafqFAiZcPzDyucfv/3pA1g88SeXUA9vxtw8WMKCHV2Lp2mYvrui/HxccGd6I/NiwI/hixQYahpQoYtUfEb+R4G/nFwN7cNR+rUL7lPsG7F+18b4KRGgd0th+eESPIz92+UFC8OqXw/g95vZu1iB08En/YiWsJJFJUL/02C7c+tosWjeJrmKsX1S57/dJ6Ixx/WkHcP1pqR+T2cndd7HGPUq3Mqhr3Htfqh55iV1FU7VBhGnbpAFTbx5EkyqObdS4flFWjdup+KUmv5djsCSRrSG99+acvh3TJ8yRBYkMiAjfOyz1l3HzGb34z5SvayhHlc44ZB/67Nsi6Yecjn91FaZRcfhh8frPjqPVXvEnljMP3SdtgDqsc4u4UpOIcEOKk0yzRvUie2g9dVk/1m3ewYWPTXRTEqvdKk8+bZrU5zsH7xM3P/G32LZpAxaPOj20bSUV/8d9Zf/9ePhD797QxJLhr07uwX3vzqWoQGInxViQSOgZV7+ogG3lmTWq+1f/lx+/X1JPpTAn9GjDK58vi1VXXnJsFy6ppi6bifxB5/wLAj8w1vRdNs0aZnZPQk2YdsugWO/An53UncmLJyV1pc/EQ2mqtPPF2iTy5JJju8T1GqpJ2QaIqurdoRntm8Vv60/nHhZXDK5ujYqL6NA8uoG0pvjF/VRX1H6DYmGBxK4c/bvpj+3WKq6KIaqaqm1IacV3Qo82cZ0jonoEnXlYB6bePCjr0mY+NKhXyPw7hsTusfA/c7ZtR6cfVHdK8blq2qBerAr02G6tWXDnaRnfWFcbLEiYXZLf2J6qV1Y+BDsLBGVykvMfJ3nuEZ3o36MNb/+qPzedfgATrh/IxceUxFU3RXU7DVb7gNeOExTsIHF+v+h7QGrzSrqosCBQ9+6qm7KsU7//vMOYfdvgvOfNpGfVTWaXdP2Qnlx+fJdYr6uHLjicdZt3pFmqMqg0aVCU1OsqTIfmDeK6LPqC57jTD27PG9OWJ6WpV1jArJGn0sCd6Hu4XjJ+VV+wJJHJXeHTbxmU1AjfonExR3ZpycRFa6u9l0s+HNu1NS9/tiwp2KXjDWpYPXcUm9QsSJhdUlFhQVzV15AMqyO6td2Ll646mt4dmsVdpb/5i+NZuWErF/9jckbrCZYkHjjvMP4SMRxLVPsOxHeh9U/wqbr/1i8qTLobf1fz/T4dOaln21gXU1P3WZAwe5w++ybf5HhA+6ahVUvDjinhyzWbk24KDF60i0jam/jCBHsu+Yu/938DItNHFTb8ULMLFCQALEDsYnbtyxJjqlnzhsVJwyRAfm7aO6prcg+uVLVONXWjoDFBFiTMbsnvTZQrkfDnjeTjhO0/PyR+g6nzEsqevmuqkVU3md3O27/qn9UduamIRA0amI91J485lepZIlEN0/5wE1bOMNXBShJmt9OjXZO8dfkscM9uSJSvnkQ3n9Er5UCQ2dgVejeZXU9OQUJEWorIOyIyz/0PHZNARIa5NPNEZFhgeh8RmS4i80XkfnFHuYjcIiLLRGSK+0v9ZBNjqkmBSLU+3vaSY7vEDa9dleHHq/lJqmYPl2tJYjgwTlW7A+Pc+zgi0hK4GTgS6AfcHAgmDwE/Brq7v+DdMvep6qHub0yO+TSmSkTyU7WUwZYi53wvw+c9W0HCVIdcg8RQ4En3+kngzJA0pwLvqOpaVV0HvAMMFpH2QFNVnaDeIDZPRSxvTK0pkJodcTOsIHHPOYfwxUi729jUjlyDRDtV9W81XQG0C0nTAVgSeL/UTevgXidO910jItNE5PGoaiwAEblCREpFpLSsrKxKH8KYKKo10/U01SYKCyR0RFzffT88lPOP7MxhEc+mNiYXaYOEiLwrIjNC/oYG07nSQL5qRx8CugKHAsuBP0YlVNVHVLWvqvZt06Zqwz4bE6W8QmukuulX7nnZVRnorVPLRtz5vYOq7RnHZs+Wtgusqp4cNU9EVopIe1Vd7qqPVoUkWwYMCLzvCIx30zsmTF/mthl7oomI/B14PV0+jakO5Tu1Rqqbzj+yM+dHPF/DmNqU66XHaMDvrTQMeDUkzVhgkIi0cNVGg4Cxrppqg4gc5Xo1XeQv7wKO73vAjBzzaUyV7KiosDudzR4t15vpRgEviMhlwJfADwBEpC/wE1W9XFXXishtgD9y2khVXeteXw08ATQE3nR/AHeLyKF41VeLgStzzKcxVbKzQqu1C6wxdV1OQUJV1wADQ6aXApcH3j8OPB6RLumJ4qp6YS75MiZXvx3ck7KN2zi8c4vQR4gas6ewYTmMCbFP8wZcNcB73sTOCitJmD2XBQlj0gg+EOhHR3VmydottZgbY2qWBQljQgQbq4Odm24/86BayI0xtcc6VhsTIvGhQsbsqSxIGBMi1ZDdxuxJLEgYE8IKD8Z4LEgYE8JihDEeCxLGhLB2CGM81rvJmBCJMeLG0w7giDw9Qc6YXYkFCWNCJJYjftx/v1rJhzG1zaqbjAlh1U3GeCxIGBPCQoQxHgsSxoQosF+GMYAFCWNC2c10xngsSBgTxmKEMYAFCWNCWYwwxmNBwpgQ1rvJGI8FCWNCFFiMMAawIGFMKGu4NsaTU5AQkZYi8o6IzHP/W0SkG+bSzBORYYHpd4jIEhHZlJC+vog8LyLzRWSiiJTkkk9jsmW1TcZ4ci1JDAfGqWp3YJx7H0dEWgI3A0cC/YCbA8HkNTct0WXAOlXtBtwH3JVjPo3JisUIYzy5BomhwJPu9ZPAmSFpTgXeUdW1qroOeAcYDKCqE1R1eZr1vggMFGtJNDXJjjZjgNyDRLvASX4F0C4kTQdgSeD9UjctldgyqloOrAdahSUUkStEpFRESsvKyrLJuzGRCuyaxBggg1FgReRdYO+QWTcG36iqiojmK2OZUtVHgEcA+vbtW+PbN7snCxHGeNIGCVU9OWqeiKwUkfaqulxE2gOrQpItAwYE3ncExqfZ7DKgE7BURIqAZsCadHk1Jl+sdtMYT67VTaMBv7fSMODVkDRjgUEi0sI1WA9y0zJd79nAe6pqpQRTYyxGGOPJNUiMAk4RkXnAye49ItJXRB4FUNW1wG3AZPc30k1DRO4WkaVAIxFZKiK3uPU+BrQSkfnAtYT0mjKmOlmMMMaT05PpVHUNMDBkeilweeD948DjIemuA64Lmb4VOCeXvBmTC6tuMsZjd1wbE8JihDEeCxLGhLAYYYzHgoQxIay6yRiPBQljQliIMMZjQcKYEHbHtTEeCxLGhLAYYYzHgoQxxphIFiSMCWElCWM8FiSMCWFPpjPGY0HCmBAF9sswBrAgYUwoK0kY47EgYUwIa5MwxmNBwpgQFiOM8ViQMCaElSSM8ViQMCaEjd1kjMeChDEhLEQY47EgYUyAX4CwkoQxHgsSxgRIwn9j9nQ5BQkRaSki74jIPPe/RUS6YS7NPBEZFph+h4gsEZFNCekvFpEyEZni/i5PXqsx1ccKEsZ4ci1JDAfGqWp3YJx7H0dEWgI3A0cC/YCbA8HkNTctzPOqeqj7ezTHfBqTFRsq3BhPrkFiKPCke/0kcGZImlOBd1R1raquA94BBgOo6gRVXZ5jHowxxlSTXINEu8BJfgXQLiRNB2BJ4P1SNy2d74vINBF5UUQ6RSUSkStEpFRESsvKyjLOuDFh/AZrK0gY40kbJETkXRGZEfI3NJhOVRXQPOXrNaBEVQ/GK3k8GZVQVR9R1b6q2rdNmzZ52rzZ01nvJmM8RekSqOrJUfNEZKWItFfV5SLSHlgVkmwZMCDwviMwPs021wTePgrcnS6fxuSD9W4yJl6u1U2jAb+30jDg1ZA0Y4FBItLCNVgPctMiuYDj+y7wRY75NCYr1nBtjCfXIDEKOEVE5gEnu/eISF8ReRRAVdcCtwGT3d9INw0RuVtElgKNRGSpiNzi1vtzEZkpIlOBnwMX55hPYzJSeTNd7ebDmLoibXVTKq5aaGDI9FLg8sD7x4HHQ9JdB1wXMv164Ppc8mZMVXjPkVCrbjLGsTuujQmyRglj4liQMCagMkZYlDAGLEgYE6rAYoQxgAUJY+LYKLDGxLMgYUyAX81kIcIYjwUJYwLUDRpQYPVNxgAWJIyJU1zo/SQaFRfWck6MqRtyuk/CmN3Nv39yDB/OLaNeoV0/GQMWJIyJs//eTdh/7ya1nQ1j6gy7XDLGGBPJgoQxxphIFiSMMcZEsiBhjDEmkgUJY4wxkSxIGGOMiWRBwhhjTCQLEsYYYyKJqtZ2HvJGRMqAL6u4eGtgdR6zky+Wr+xYvrJXV/Nm+cpOLvnaV1XbhM3YrYJELkSkVFX71nY+Elm+smP5yl5dzZvlKzvVlS+rbjLGGBPJgoQxxphIFiQqPVLbGYhg+cqO5St7dTVvlq/sVEu+rE3CGGNMJCtJGGOMiWRBwhhjTKQ9MkiIyDkiMlNEKkSkb8K860VkvojMEZFTA9MHu2nzRWR4DeTxeRGZ4v4Wi8gUN71ERLYE5v2tuvOSkK9bRGRZYPunBeaF7rsaytcfRGS2iEwTkVdEpLmbXqv7y+WhRo+dFPnoJCLvi8gsd/z/wk2P/E5rMG+LRWS6236pm9ZSRN4RkXnuf4saztP+gX0yRUQ2iMgva2N/icjjIrJKRGYEpoXuH/Hc7463aSJyeE4bV9U97g84ANgfGA/0DUzvBUwF6gNdgAVAoftbAOwHFLs0vWowv38ERrjXJcCMWtx3twC/Dpkeuu9qMF+DgCL3+i7grjqyv2r12EnIS3vgcPe6CTDXfW+h32kN520x0Dph2t3AcPd6uP+d1uL3uALYtzb2F9AfODx4LEftH+A04E1AgKOAiblse48sSajqF6o6J2TWUOA5Vd2mqouA+UA/9zdfVReq6nbgOZe22omIAD8A/lUT28tB1L6rEar6tqqWu7cTgI41te00au3YSaSqy1X1M/d6I/AF0KE28pKhocCT7vWTwJm1mJeBwAJVreqIDjlR1Q+BtQmTo/bPUOAp9UwAmotI+6pue48MEil0AJYE3i9106Km14TjgZWqOi8wrYuIfC4iH4jI8TWUj6BrXDH28UAVQG3uo0SX4l1J+Wpzf9Wl/RIjIiXAYcBENynsO61JCrwtIp+KyBVuWjtVXe5erwDa1UK+fOcSf6FW2/sLovdPXo+53TZIiMi7IjIj5K9WruLCZJjH84g/OJcDnVX1MOBa4FkRaVqD+XoI6Aoc6vLyx3xuO4d8+WluBMqBZ9ykat9fuxoR2Qt4Cfilqm6gFr/TgONU9XBgCPBTEekfnKlePUqt9NcXkWLgu8C/3aS6sL/iVOf+KaqOldYFqnpyFRZbBnQKvO/oppFiepWly6OIFAFnAX0Cy2wDtrnXn4rIAqAHUJprfjLNVyB/fwded29T7bsayZeIXAx8BxjofjQ1sr/SqPb9kg0RqYcXIJ5R1ZcBVHVlYH7wO60xqrrM/V8lIq/gVdOtFJH2qrrcVZesqul8OUOAz/z9VBf2lxO1f/J6zO22JYkqGg2cKyL1RaQL0B2YBEwGuotIF3dVca5LW91OBmar6lJ/goi0EZFC93o/l8eFNZAXf/vBus3vAX5vi6h9V1P5GgxcB3xXVTcHptfq/qL2jp0krn3rMeALVb03MD3qO62pfDUWkSb+a7xOCDPw9tMwl2wY8GpN5isgrjRf2/srIGr/jAYucr2cjgLWB6qlsleTLfR15Q/vi12Kd4W5EhgbmHcjXm+UOcCQwPTT8HqDLABurKF8PgH8JGHa94GZwBTgM+CMGt53TwPTgWnuYGyfbt/VUL7m49XDTnF/f6sL+6u2jp2IfByHVyUxLbCfTkv1ndZQ36hTrAAAAilJREFUvvbD6/U11X1XN7rprYBxwDzgXaBlLeyzxsAaoFlgWo3vL7wgtRzY4c5dl0XtH7xeTQ+64206gR6cVfmzYTmMMcZEsuomY4wxkSxIGGOMiWRBwhhjTCQLEsYYYyJZkDDGGBPJgoQxxphIFiSMiSAizUXk6sD7fUTkxWrYjj/09MgUabqKNyz1pnxv35hU7D4JYyK4QfBeV9Xe1bydW4BNqnpPBmk3qepe1ZkfY4KsJGFMtFGAfwX/B/EeYDQDvHGiROQ/4j3sZbGIXCMi17rRZieISEuXrquIvOVGN/2viPRMt1EROUEqH2jzuT9khTG1Ybcd4M+YPBgO9FbVQyFWsgjqjTfcdgO8YUF+q6qHich9wEXAn4BH8IZWmSciRwJ/BU5Ks91fAz9V1Y/diK1b8/R5jMmaBQljqu599R7es1FE1gOvuenTgYPdCf4Y4N/e2HqA9+S+dD4G7hWRZ4CXNTDAozE1zYKEMVW3LfC6IvC+Au+3VQB845dEMqWqo0TkDbzB9z4WkVNVdXY+MmxMtqxNwphoG/GeBV0l6j3QZ5GInAOxB9Qfkm45EemqqtNV9S68ocbTtmMYU10sSBgTQVXX4F3JzxCRP1RxNRcAl4mIPwx2Jk9G/KXb5jS8oaHfTLeAMdXFusAaU8usC6ypy6wkYUzt2wRckcnNdHgPyTKmxlhJwhhjTCQrSRhjjIlkQcIYY0wkCxLGGGMiWZAwxhgT6f8BRpY7L25aJO4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# find the segments of good data for both source and receiver\n",
    "bb=np.intersect1d(sou_ind,rec_ind)\n",
    "if len(bb)==0:raise ValueError('Abort! no good data in overlap')\n",
    "\n",
    "# do cross correlation\n",
    "corr_day,t_corr,n_corr = noise_module.correlate(source_white[bb,:Nfft2],receiver_white[bb,:Nfft2],fc_para,Nfft,dataS_t)\n",
    "\n",
    "# plot the waveform\n",
    "print(len(corr_day))\n",
    "tvec = np.arange(-maxlag,maxlag+dt,dt)\n",
    "plt.figure()\n",
    "plt.plot(tvec,corr_day)\n",
    "plt.xlabel('time [s]')\n",
    "plt.title('cross correlation function between AYHM and ENZM')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Save cross correlation data into ASDF file\n",
    "\n",
    "Though we only have one station pair, we can try to save it into the ASDF file. We save the cross correlation data into the auxiliary structure of ASDF, which has two dimentions (data_type and path). In this example, we use the station and network name of the source and receiver station to define the $data\\_type$ and use the channel names to define the $path$. The two tags are chose because the two-dimention variable are enough to define any cross component of the cross correlation functions for any station pairs. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "cc_h5 = 'cc_example.h5'         \n",
    "with pyasdf.ASDFDataSet(cc_h5,mpi=False,mode='w') as ccf_ds:\n",
    "    # location info \n",
    "    coor = {'lonS':slon,'latS':slat,'lonR':rlon,'latR':rlat}\n",
    "    # cross component\n",
    "    comp = tr_source[0].stats.channel[-1]+tr_receiver[0].stats.channel[-1]\n",
    "    # parameters to be saved into ASDF\n",
    "    parameters = noise_module.cc_parameters(fc_para,coor,t_corr,n_corr,comp)\n",
    "\n",
    "    # data_type name as source-receiver pair\n",
    "    data_type = tr_source[0].stats.network+'.'+tr_source[0].stats.station+'_'+tr_receiver[0].stats.network+'.'+tr_receiver[0].stats.station\n",
    "    # path name as cross component\n",
    "    path = comp\n",
    "    # command to save data and parameters into asdf structure\n",
    "    ccf_ds.add_auxiliary_data(data=corr_day, data_type=data_type, path=path, parameters=parameters)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Read the ASDF file\n",
    "\n",
    "Finally, we want to read the cross correlation function we just saved. To retrive the data, we simply need the two tags we just created for the auxiliary structure in ASDF, which are $data\\_type$ and $path$. Note that we do not necessarily need to know the two parameters beforehand, because we can simply get the two parameters from reading the file. You will see how we do it from the codes below.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['E.AYHM_E.ENZM'] ['UU']\n",
      "[-0.00122514  0.00055406  0.00109561 ..., -0.00040456 -0.00018136\n",
      "  0.00081051] {'azi': 185.50987, 'baz': 5.5054541, 'cc_method': 'xcorr', 'comp': 'UU', 'dist': 7.1563296, 'dt': 0.050000000000000003, 'latR': 35.60844, 'latS': 35.672642, 'lonR': 139.70786, 'lonS': 139.71544, 'maxlag': 100, 'ngood': 156, 'substack': False, 'time': 1292457600.0}\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAD6CAYAAABUHLtmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd5wU9f348df77jiadI4i7agiYoMTURQLiGAJ9lhisAVjiUn8JQYbGixRY/QbE2MkStQkBo0loCIgKBpjKIfSpRxF4aQcRYrUu3v//tiZvdndmS23ewV4Px+PfdzuZz4787nZ3XnPp8xnRFUxxhhj/GTVdAGMMcbUXhYkjDHGBLIgYYwxJpAFCWOMMYEsSBhjjAlkQcIYY0ygjAQJERkqIstEpEhERvksrysirznLZ4lIvpN+jojMFZGFzt+zPe/p66QXicgzIiKZKKsxxpjkSbrXSYhINrAcOAdYB8wBrlLVJZ48twLHqeqPReRK4GJV/b6InAhsVNVvRKQ3MEVV2znvmQ3cAcwCJgHPqOr78crSsmVLzc/PT+v/McaYw83cuXM3q2qe37KcDKy/H1CkqqsARGQ8MBxY4skzHHjQef4G8EcREVX9wpNnMVBfROoCzYHGqjrTWecrwEVA3CCRn59PYWFh+v+RMcYcRkTkq6BlmWhuages9bxe56T55lHVUmA70CIqz6XA56q6z8m/LsE6jTHGVLFM1CTSJiLHAI8DQyrx3pHASICOHTtmuGTGGHN4y0RNohjo4Hnd3knzzSMiOUATYIvzuj3wNvBDVV3pyd8+wToBUNWxqlqgqgV5eb5NasYYYyopE0FiDtBdRDqLSC5wJTAxKs9EYITz/DLgQ1VVEWkKvAeMUtX/uplVdT2wQ0T6O6OafghMyEBZjTHGpCDtIOH0MdwOTAG+BF5X1cUiMkZEvudkexFoISJFwJ2AO0z2dqAbMFpE5jmPVs6yW4EXgCJgJQk6rY0xxmRe2kNga5OCggK10U3GGJMaEZmrqgV+y+yKa2OMMYEsSBgTZdt3+3lvwfqaLoYxtYIFCWOi3PRKIbe9+jlbdu2r6aIYU+MsSBgTZfnGnQBk2XRhxliQMCZaefmhM5jDmHRZkDAmgIUKYyxIGBPDnZW+/BAaHm5MZVmQMCaKe+2QxQhjLEgYE0itwckYCxLGRHObm6wmYYwFCWMCWZAwxoKEMTHcqyOs49oYCxLGxNCov8YczixIGBPgUJoh2ZjKsiBhTBS3uclihDEWJIwJZEHCGAsSxsRyqhLWcW2MBQljYmnEH2MOaxYkjAlgHdfGZChIiMhQEVkmIkUiMspneV0Rec1ZPktE8p30FiLykYjsEpE/Rr1nhrPOec6jVSbKakxC4eammi2GMbVBTrorEJFs4FngHGAdMEdEJqrqEk+2G4FtqtpNRK4EHge+D+wF7gd6O49o16hqYbplNKZyLEoYk4maRD+gSFVXqep+YDwwPCrPcOBl5/kbwCAREVX9TlU/JRQsjKkVKq64rtFiGFMrZCJItAPWel6vc9J886hqKbAdaJHEuv/qNDXdL+J/L0kRGSkihSJSWFJSknrpjYkSvuLagoQxtbrj+hpVPRY43Xlc65dJVceqaoGqFuTl5VVrAc2hzaYKNyYzQaIY6OB53d5J880jIjlAE2BLvJWqarHzdyfwKqFmLWOqXLi5qbxGi2FMrZCJIDEH6C4inUUkF7gSmBiVZyIwwnl+GfChxhlfKCI5ItLSeV4HuABYlIGyGpM0q0kYk4HRTapaKiK3A1OAbGCcqi4WkTFAoapOBF4E/iYiRcBWQoEEABFZAzQGckXkImAI8BUwxQkQ2cA04C/pltWYVFifhDEZCBIAqjoJmBSVNtrzfC9wecB78wNW2zcTZTOmsixIGFO7O66NqVHW3GSMBQljYrijre06CWMsSBgTyOZuMsaChDGBLEQYY0HCmEBWkzDGgoQxgSxGGGNBwpgYYlOFGxNmQcKYANbcZIwFCWMCWYgwxoKEMYHKrSZhjAUJYwJZjDDGgoQxQSxGGGNBwpgYFbcvtTBhjAUJYwJYjDDGgoQxgawmYYwFCWMCWYgwxoKEMcEsShhjQcKYaBX3k7AoYUxGgoSIDBWRZSJSJCKjfJbXFZHXnOWzRCTfSW8hIh+JyC4R+WPUe/qKyELnPc+I+8s1pppYjDAmA0FCRLKBZ4FhQC/gKhHpFZXtRmCbqnYDngYed9L3AvcDv/BZ9XPAj4DuzmNoumU1JhVWkzAmMzWJfkCRqq5S1f3AeGB4VJ7hwMvO8zeAQSIiqvqdqn5KKFiEiUhboLGqztTQLGuvABdloKzGJM1ChDGZCRLtgLWe1+ucNN88qloKbAdaJFjnugTrNKZKWUXCmEOg41pERopIoYgUlpSU1HRxzCHA7fyyqcKNyUyQKAY6eF63d9J884hIDtAE2JJgne0TrBMAVR2rqgWqWpCXl5di0Y2JpVF/jTmcZSJIzAG6i0hnEckFrgQmRuWZCIxwnl8GfKhxTtNUdT2wQ0T6O6OafghMyEBZjUmadVwbAznprkBVS0XkdmAKkA2MU9XFIjIGKFTVicCLwN9EpAjYSiiQACAia4DGQK6IXAQMUdUlwK3AS0B94H3nYUyVq2huqtFiGFMrpB0kAFR1EjApKm205/le4PKA9+YHpBcCvTNRPmMqw2KEMYdAx7UxmeZetmkd18ZYkDAmhhsbLEYYY0HCmEBqDU7GWJAwJprb3FReXrPlMKY2sCBhTACrRxhjQcIYHzZVuDEuCxLGxNCIP8YczixIGBPAOq6NsSBhjA+3uamGi2FMLWBBwpgA1iVhjAUJY2KEh8BalDAmM3M3GXMomLJ4Ay9+urriiuuaLY4xtYIFCWMcN/9tLgDNGtQJJVhNwhhrbjImmtthbR3XxliQMCZGuRMdbBZYYyxIGBPD7bC2EGGMBQljYpQ5QcKam4yxIGFMDHf2V2tuMsaChDEx3JqExQhjMhQkRGSoiCwTkSIRGeWzvK6IvOYsnyUi+Z5ldzvpy0TkXE/6GhFZKCLzRKQwE+U0Jhllbsd1VK9E8bd7UFUefncJRZt21kTRjKl2aQcJEckGngWGAb2Aq0SkV1S2G4FtqtoNeBp43HlvL+BK4BhgKPAnZ32us1T1BFUtSLecxqTKW5NYVLydAY99yOOTl/HCp6u57q9zaq5gxlSjTNQk+gFFqrpKVfcD44HhUXmGAy87z98ABomIOOnjVXWfqq4Gipz1GVPjvB3X33y7B4AZyzYB1hRlDh+ZCBLtgLWe1+ucNN88qloKbAdaJHivAlNFZK6IjAzauIiMFJFCESksKSlJ6x8xxuvb3ft5cOJi9peWI86ETtUdHHbvL63eDRoTpTZ3XJ+mqn0INWPdJiID/TKp6lhVLVDVgry8vOotoTmkPf/JKl76bA0T5hU7k4dX7z0mFhVvp9foKby/cH21bdOYaJkIEsVAB8/r9k6abx4RyQGaAFvivVdV3b+bgLexZiiThLVbdzN/7bcZXWe5anhm2OqsSSwq3g7AR04TlzE1IRNBYg7QXUQ6i0guoY7oiVF5JgIjnOeXAR9qaBD6ROBKZ/RTZ6A7MFtEGopIIwARaQgMARZloKzmEHf6Ex8x/Nn/Zny94SCR8TUHy8kO/TxLfa7q+8esr3hvgdUwTNVLexZYVS0VkduBKUA2ME5VF4vIGKBQVScCLwJ/E5EiYCuhQIKT73VgCVAK3KaqZSLSGnjbaQfOAV5V1cnpltWYylAl3CdRtGlX3Lzl5crOvaU0cWeSTYMTI8JDcr3ufTt0znT+ceenvR1j4snIVOGqOgmYFJU22vN8L3B5wHsfAR6JSlsFHJ+JshmTCZI4CwD/N205z3xYxNKHhlKvTnbiN8SR5QQmvyDhtX33AerlZlE3J73tGeOnNndcG1MrKBU1CZcERI0XP10NwHf70h+V5AaJRE1cx4+ZylVjZ6a9vUxTVe7790IWrMtsH5GpXhYkzGHj1n/MZcK86DEVianG1iTWbdvjmzc7K5Rzb2loAqjXC9eSP+o9djlBY39pORPnf5NwXqhNO/fy7Z4DzvYT94R8/nXtOxDv2FPK32d+zQ9emFXTRTFpsCBhalx5ubJ683dVsu6vtnwXDgyTFm7gp+PnVWo9QTWHaLlOk8/eA2UAPP/xSgDWOxfjPTN9BXf88ws+WLIx7nr6PTKd+/99aIzVsOsOD24WJEyN+9OMIs56cgbLNiQ/H9Lqzd+xZde+hPmG/f4/lQ4MLkVZuj65suVmOzUJJ0hIVJPRN9tDwWL7ngP8bPwXFDw8LfH2D9ajbMXFJeYgZkHC1LjCr7YBUPzt7qTfc9aTMzj9iY8S5tu9vyzu8tKy8vABPciCtdt5ZNKXSZXLDQqlZaEjo3uc/OLrbRR/W9FEdaBM+fe8b9i8ax+PT17Kxh17A9fpDRKbd+1j9ISDo4aRbO3L1G4WJEyNC59wpnjGuXt/Gb9+Z3Fa277qLzPpeX/80dWvFa6Nu9yPO924e6D81ZsLGfDYh+Gz6gcmVhzon5uxkl/8a37gurxXeT/87hJe+d9XKZenJqVakZi5agvzMnxB5KHm4+Ul4XnEqpoFiRqwfc8Blm7YUdPF8LW/tDzhkMtMS2depL/+d01a256zZlta7w8ya9VWACSqy9v9Fw+URf6z+0vLmbNma3giwSAH493yUr1505VjZ3JRFVwQGaRo0042xanJVac9+8tYWRL/WhyAEeNmV9tMxBYkkjB79VZO+c108ke9R2lZedrru3LsTIb+339Ses/ib7Yz6s0FlFfiKLFnfxnbvtvvu6ysXNlfWvE/9bjvfa7+S+LhlCU797FuW/LNQ/Hs3OuM4glYvvdAGd3umRQ4MsnvILQ14P8Nsnxj/D6HGwZ0jkn74uttEQeX7/aVst0ZkfT45KVAbJPL218Ej666/M//46wnZ8Ske/+9rIOoCacyQX9VEgdIgM+KNvPcjJXOdpQ7X5vHzFVbUt8gMPipT+j/m+mVem861m7dHfOdvu3Vzxn0u4+TPs7sPVBWqWNCKixIJOGK5//H+u2hg8G+0vSDxJfrU69F3PDSHMbPWcumnYk7a6Nd/Kf/cuJDH/gu++G4WfS4733KypUdzsF61uqt/KtwbfjgHW399j2c9Mg0Tns8cZ9AMtyz+aAzzpKd+ygtV56YvMx3+Y69pTw1dRl3vVHRZLM/6nPyrnv15u8456nIH+KQpz+JW8ZWjevGpF38p88Y/NTH4ddn/HZGeKhrqvY4/SJ+36/yiCCROEqoKl9viQzgO/YeIH/Ue7wz/5tKlS8Zn63cTP6o91jhBlyN+JOUs3/3ceJMwNUvzAoH4gNlyltfFKc11LYmamgXPfvfmEEVHy8PzWRdlmSE7Xn/ZO6r4j4qCxJRVJXHJy8NHGmTyndp74EyDmSg5gGwcUdscHhyyjImJTFD6NKo/+Wn47/gppdDVdX/FoXOvkZPWMRxD04N5/nlGwu4L2AI5vKNic/2Jswr5qF3lyTM51XZ3+nTH4Sucn69cF04rbS8PKK/wnsQOOvJGazYtIu731qY9DbKA360O/ZWBIXNPqOtoi/CC7Jg3fbw88mLoj/Tim1nRVUl/ALr85+sYuBvPyJ/1Hs89O4S9peW89Xm3c6ylRwoK+f/vT6f8bO/Zu3W3SzbsDPpM/h43nXmkpq5eqtT6tRuA7vkm8o1wbqfjTeA9ntkGi/9d3X4dWlZeaVqvjOWbfL5PIJNmFfMtS/GBqudew/w1ueh76eqUl6ubHFqu97PMHxXxKh9tqh4e8T/4/XqrK+TLl9lZGRajkPJjj2lPDdjJeNnf80Xo4fEnJEm8uX6HWzbvZ9Tu7ak5/2T6dmmEZN/5jvLORAKJDlZEp7MLZH+v5nOuz85jd7tmvDHj4oAWPNYcvP3bNqxl1aN6zFhXuzZ5L99mkFKAmot8Q574z5dzTff7uEF58rj+84/mhc/XU2/zs3JzcmiZ5vGge+9+W9zuapfR35zybHx/5EofqOTyssj+ysWFm+PyfPVluQPGvGq9Lv2lVK4ZmtMev6o95Jev9f/TVsRuCy6ucmv/8i96tt93qP1EfRq2yScNmfNVt78fB1vfr4u4n3Jfo+iqSpz1mwLH+zcMqbS3KSqnPdMak2wACNfKeQnZ3cHQk1767fv4dTHPkQVHnxnCdc5zYRPTFnG2E9WMfueQbRqXC/hej9ZXsIPx80Ov05237g1g8+KNnNqt5bh9HveXsQ787+he6tGzF6zNeIEqlwhO+pzjd53F/zhUwAG92pNeWbOO5NmQcKxacdeGtWrmJTN/fFFn318u3s/OVkSnpdn23f7Kf52D73bhX6Ew34f+qK7XyrvWfymnXt54T8VP2BVpef9kzn3mNY8f23sHVq37z5Abk4W9XMj5+SZOP+b8PZS0e/R6Uy78wzfZck0Y7iCspbs3MeYqNrD4m928PB7FcNH3f3y3b5SGtbNYV9p5AH+n7O/5rxj23B69zwmzCumRcO6dGrRIO52/c7Wt+2O7JOY+1VsB3WyVXqAVXEu9vv5a/MSXhyXiuhiua/nfrWVnXsjm7P8Yld0cPdWZpdv3EXdnMo3IOSPeo+bTuvMfRf04p63F0acxTZ1JjV0O+vdou05UMb7C9fz5YadTFq4nnd/clr49/PZys00rleHHq0bxWzr8clLuXlgF5o2yGXGsk08PW0Fb91yaviqdoCpSzYy1dn3IjBjWUnE/nt9zlouL2jP2E9WAbB51/6kgoQ3QECo/6BD8wYJ3+e6+oVZNMjNZsmYoQDhvqsL//gpR0X9r2XlGvE/QfB9SzLVxJsKa25y9Ht0OkePnsz0paEvnPsRRX94pz3+UcTIi0ue+ywc5T9cWnGg8HacfrevlKmLN9DvkenhLytU/MCnLPY/wBw/ZipHj54cDjwuv2NlebkyK4mOuxFRX37XTp+29KJNu3zb2Ddsjx0Joqq+NY81WyIPrrf+Yy6zVm3hmAem8PHyEh5/P7af4doXZ1Oycx8/HT+PH7w4K3zAX7dtD/mj3mP77si+Er/gET1duF/TVyp9wG99HtzhnMxolFQsi+pEV0Kf76XP/Y/3F22IWPb76cvDz7fvOeDb/OT9Cu8vLednryV/ceGe/WV8uHQjqhpuTnNridHNHN96PpcpizfQx9MPdss/PueZ6Sso2rSL6V9uose97zNv7bdc/ZdZXPCHT32vE3luxkpOGBNax6/eXMD8td/GvZ5EkJgAe9ebC5i0sGKflZVrRZ9JCrzX5KhqUhd+7t5fxurN3/H3mV8xa3VFTTP68/Ub0VabRrFZTSLKna+HOj937i3l6y27Y4IEVNQONmzfGzGdxA0vFYafe38gxzwwxXdbQe3c0aI7uj9eXhLTwfnCp6t4dNJSrihoz73n9WL3gVJaNKxLbtRZY3GCIZZem3bu4/q/zmbTzn00rV+HCbefxqSF6/nlGwvCeXbvLyVLhNcL1zJ6Quw1C7e/+kXE60kLN4R/tJ8Vbebrrf5n6Od7mh7+7Ext4Xp3YWRz2a69qXUWD+yRxyfLS/j+SR3CF/LVZqoaWOt59qOKfXP8r6fyxKXHxeSZMO8bjm5b0cy3dmv878CBsnLeW7Ce4SccydGj/a8hmf5lcM1JBH4fp8ns06IS9peVR5xsJbow0q15xGv+3XOgjGedJlgvd8QZhM7kASb/7PTAps+gkUXrt++hbZP6PPXBcv7wYRGv33wK/To3Z9SbCxg/Z61vk5TfaLVoZz45gzWPnc/arRWtFu6xYde+0pS/35lmQSKOgb/9iJ+c3c132dINO1IexhrN256890BZ0lNLL92wM6Yz2u30e71wHV9v3c3MVVu55MR2lWqW8nJHHn0FrNi4k1v/8XnE8l6jQwHw9O4to9+a0POeWlU07ygu75kghIYke01MccTO107tJjqAVtaqkqqZd8r13b4yhv8xuesG3lkQuy/+t2oL16Qw8uf5j1fy5NTl4RFXfm58uTBwWaIBAaleh3Pyo9PCAzfemf8NPxnUPTCv30lQdJMmwPpv91K8bY/v//Hk1OUxaQCn/OZD5o0+hz98GApEL/xnFb95/0u+cCZX9PYFpWp/aXlEh/dNLxfy5GXHM/C3yTUvjf1kJSMHdq309uOx5ibgo6XBVy66X4ho0QHiXZ8fZyLes+We909m6uIN7NpXyuZd+7jsuc+SXs+cNVs54PnhFW0KHbTe+qI4po/AK9Wrlc+JM0z0Pys2p7SudKQ7YmyN02Gd7pxO1WX2mq0sSXLY9OZd/teHJDM0122qcg+SqYz+SoV3FFoyvCP7fveB/wE8nl+/4/8biA4Q5eXKjGWbmPtV7CAE1zPTK44HU5dsDAcI8G/STFaP+94Pfy8hdCK0cnPyzZiPTlrK+u3JtxKkwmoSEB4llI7oZpVkrIw6Ax35t7mV2vZTU5fzP09/hN9QTD/pXq1cU6JrFqZCZa7BcXW9ZxK/Ht47g6WpGu5FdOnwO9Hocs8kn5yRxgUMQ60KZWWp1bhe+mwNdw87OuPlkFQvma/NCgoKtLAwuBoc5LLnPjso2qaNMYePOtkSM31LPCflN+NfPz61UtsSkbmqGjvEkgw1N4nIUBFZJiJFIjLKZ3ldEXnNWT5LRPI9y+520peJyLnJrtMYYw5lqQQISG0Ye0rrTXcFIpINPAsMA3oBV4lIr6hsNwLbVLUb8DTwuPPeXsCVwDHAUOBPIpKd5Doz5tCpSxljDlc50VfkZUgmahL9gCJVXaWq+4HxwPCoPMOBl53nbwCDJHQF1HBgvKruU9XVQJGzvmTWmTGltWlQsjHGVEKtrUkA7QDvhPvrnDTfPKpaCmwHWsR5bzLrzJhMzOxqjDE1ye+arkw46IfAishIESkUkcKSkpJKrSNTk/AZY0xNya7FNYlioIPndXsnzTePiOQATYAtcd6bzDoBUNWxqlqgqgV5eXmV+gdKU+wgMsaY2qZh3aq5oiETQWIO0F1EOotILqGO6IlReSYCI5znlwEfamjs7UTgSmf0U2egOzA7yXVmjPVJmIPJL889qqaLYGqh287ynx0iXWkHCaeP4XZgCvAl8LqqLhaRMSLyPSfbi0ALESkC7gRGOe9dDLwOLAEmA7epalnQOtMta5Ccg+l2X+aw175Zfa7t36mmi2Fqmcb1a29NAlWdpKo9VLWrqj7ipI1W1YnO872qermqdlPVfqq6yvPeR5z3HaWq78dbZ1W54qQOiTNVo+Papzffkql6r950co2tPztLqqyT0hy8rOO6ClVVh0+qLu/bHoChvdukva6+nZqlvY5Dxdhr+2Z8nanci6IyvDesiZYtEr7Hhp9E94s495jWMWkdU7hXQnW45MQqG8xYa8wfPaTS7+2a1zAmLSerag7nFiRI7wd/3rGRB/T2zepXel3NG+YCFTdtqaybz+jCm7dU7vL8ZNWWJrq8RpH3nl7063Nj8rj7NZNqshsrK0sYcUp+4PIBcQIM+H+/rh+QX+k701WFTM3QW5X+kWZtskmDOokzAX06No1J86s1WE2iCl14/JE8cGHsBd2tog5AAO/cflr4eevGdTnmyMimod9feaLvNnq3C81d363VEUBscIGKKZSjP+vBR7eKU/qQXw3tGX5+1lGJ8ydjYA//0WJPXXG878E4kTrZwsIHK3/25GfstX1p4Llzn9/PpCp+PF1aNkzrDm/pyBYhK0sCa4vZWcLPB/dg1LCevsu/2x87I+z5x7UF4OUb+nF8h9iDEsCd5/SISVv9m/PilrVH6yPiLg+S7LT5uVG3/W3WoA5fjhnKJSe2C98pzxV9RzjXc9f0qVQZvcH4Buc2qZk2fmR/hp8QW6vq2DyyJpGbk0XD3OT2WaosSADtmtbnmpM7cWrXFvygf0cAzu7Zik9/dXZEvmG923Csp78gSyTmTmDtmkbWJO44uxsPXNiLibedxprHzudP1/Th9O4t6dMx8gf+2sj+4VFWdaK++Jf0aR9+/s8f9ff9H7zHwf5dWsT7d30VdGpGbk5WxNnkKzf08817SZ/2MT/OVY+ex6ldW/D9guD+nZ8O6h5xi1iv1o1jA/KIUyo6Z/963Um+72tcv05Ec6HfVafZWcIlfdpxWreW4c83XVlZws8Gxx40E6nsQdPLDXpBoS8nS/jp4O5cPyDfd3n0tO5rHjufVo1Ct/Q8o0cer430/47d4XMfB79bx3qN83xuyYzKmnvfYB66qDdHtfE/oHvNG30Os+8dxMe/PDNcthGn5lM/N5unvn8C80YP4Z7zKgJlbk4WDw0/JmY9p6VwL5TPRp3tmz7a5yTz8r7tueXMrpwWp2bnF5CH9KpoDuzfpQVXn9yRW84M3SviiLo5PHdNHx67NPI+8MsfHkZOtjU3VancnCxe/VH/8AH57J6twgfNZ6+OPNN46frQF/+GAZ3JCrg3bU6W8PrNp3DnkKO43pOvR+tG/O3GkyMCwR2DunNylxZ0dw4gXfIa8lPnS//4pcdS3zmrapibzSldW/DF/efElN/vt/qfu84KV9vvO79iCuEJtw0IP+/QPBTUXhhRwPKHhwXun+iDeFaWcKKnGpyVJbz6o/48ftlxzB89hFd/FFsVdw8ofndPm3XPYP5wVWQtLL9lQ9685VTGDD+Gs3q24q1bY5vQBMj2zFnjtx+yRHjqihP4+00ns/dAZi6czJLg+xDHk52BdmP3uxQ0DYMbROpUclvxal6z7xkUfu42Ob55yynhtNdvPiUi/xGesfu3ndWNk/Lj95W1OKIu1/bvxAVOzSZI0wZ1aNogl6YNcunUIrZ93hU9vP2cXpE1+Kv6daBRvTq+LQmDj24d/q27WkfdH/uso/Ko43z/OjSvT99OzbjguLY0qpfDr4b15FdDe/L3OM1S0bs6S2DsDyMnY62TncXtzvDWclWGHds2Yr/eHVBjzBS7n0SUPh2b8ckvzwofPP2ceVSr8Bn37v2lfPPtHv4+M3S/X/fg36dTM/p1bh64Du9dv9zvydX9OnJcu6Yc274JZ/TI45Yzu1KvTja795fSL785v3bOgpo1zOX2s7pF3AdDEC44ri1n96xoaurQvAEXn9CO1wrXcnnfDpzYsRnrtu2OOHt57pq+bA5Ufc4AABXWSURBVN61j6YNgtvt2zWtz4TbB1Dw8LSI9LdvHUD+qPdi8jdpUIdTu7ake6sjWLEp9sYpBQEHiguPP5JBR7cK3+3u8oIOHFE3J9ysEl37gtDkjN7+kURjEK7q14E35qZ20xs/oVpkbHrPNo148vLjw/c9j1aegc4Mt+YU9L9e7HT6ZmUJ9553NBPmF7OouOI+E13zGsbcy8Rv/RCqYXqn0W/lOUi+d8fpAPTt1JxPfnkW7ZrVD6zdNKoXe6iZdudAPlm+mcsK2nPcg1Oj8tdh7LV9A++xEr3vB/VsxTPTV8Q0tZ57TBuemFxxH/U2Terx2aiz+WDJRk7p2oLuTvPv9QM6s3Dddt76oph//fgUuuUdQbOGuezZX/E7vbZ/p/BB3a2R/vX6itr2f+7yr2V4Hd22ccQ9P9xA//rNp3DF8/+LaDbO9wxOcAO327TqDeQjB3ZJuN10WE3CR8cWDXyr0X4/yga5OTx8UUXVr+URdXn+2r48/4P4I2qu9Ay7dZsFRCTcnCUi4XbZBrk5vP7jUyLuUxx9FisCf7y6T0TTFMBDF/Vmxi/OpEmDOvTt1CymfTOvUV3OjNOHUfTIMP5z11m0PCK2OSiRD+48gzWPnc+Pz4i8raL3OHlq1xb8YkhFs02D3JyIPoYg3h/Qby6pqJlU1SRn0UT8D/h3DOpO73ZNYjqBzzoq1L/jDpI45sjGXNKnciN4suPUJC7r255BR1c0V/xoYBd6tAo13YwZfgwLHhwScfvQj35xZsw63JrKVf068s+ApicgokmoY4sGZGeF+koe9Dkrd0vqnvUPPaYN3Vo14obTOtM4oAlyyDFtYkZiuTXs6NvlHt+hKWseOz+m+aZr3hFMvD1Uc3Z/M0c2rc+IU/Pp0bpRxO/80UuOZcrPBnJSfnOaOYMd6udmh/sybjmzKyLC8oeHMeZ7qd2c6ZGLe9Ot1RE8eXlkLdo91qsqax47n5ud38rsewbxrhOEIdRHc895PRk/MlRT8wbyRE1+6bKaRBLci1TcdttEzj0m8RDWpg1y0xpNMqx3W579aCWndm3BZyu3BObLzckiv2VwdTy6+hxKqxu+ZWQm2zkrvsuhH2vPNo141aePpVfbxhR+tS3u0GT3AKkK53jacL0Hzh/078jfZ34dtxP0iUuP4643F6TwX1Rs54qTOsTcTjNo1NepXVvy0bKS8AG6deN6NK3vX3t7747Q4IiZdw9iy3f7OP+ZyFqJ27zh15r02CXHxqS5IaFhbg6N69UJB6pPfnkWHQOG0i57eCh1srLIyhI+G3V2RPPozQO7MO3Ljb7vA7huQGcGHd2a/WXl4TN+N/C46zm9R+RB/q/Xn0S+T7PR89cWsH33AY4fE6pp/PycHlzSp53v9zZImyahvOf2iv+7rFcn27cv5K1bT+WLr7/lSKe/sTIjr645uRPXnNyJHXsPRKS7B/jo841WPv+f9x7W0c3cVcmCRBJO69aSp79/PMN6x28nrU7u2eo787/hs5Vb6HVk48Rv8ph596DA+x5/9IszU77hSTznHtOaP3+8kjN7hGosXfOOYOTALlxzsn8n8osjTuLLDTuo71OjCDd/BPxGvL+dnw7qQf8uLXzHlLvOPKpy831lidDSp7O9T8CIIzd2uUFCCLUv+3FHzLVpUs938kn3ZMWvJpFMUL9hQGd+/c4SWjYKbmKsm1Ox74+MGoxx93lHc/d58W+T2cG57mKLcyvdiqCuEa9d8UbkRQ8VjdcH4adVo3rMf2AIjSo5t1HDujkpdW7H49aa3FGO3ppEqob1bsPlBe0TZ0yTBYkkiAgXnxj/w3jgwl78e9431VSiChcefyR9OzWL+SEn4p5d+WmQ6/+1ePcnp9HiiMgDy0UnHJkwQJ3YsVlErUlEuCfOQaZJgzqBI7ReubEf23Yf4NoXZzkp0c1uFQefvEZ1ueC4IyOWR/8WWzWux5rHzvftW4nH/XHfPLALz38SmkAgumb488E9eHracnKyJHxQDAeJqJFxdXOy2FeaXKe6e/Z/0+ldYkYq+TmjRx5vf1Ecbq68fkBnrq+iIZvR3Enn3BMCNzBW91U2Teond01CdVjw4JDw6MCfnN2dOWtmxwylT8ZzCZq0M8X6JDLk+gGdI0YNVadUA0Rl9W7XhLZNIrf1f1eeGFENrmoNcnNo1zS4g7S6uNX9eGfUbodidpaEzxzdq+kHdGsR0cQQ1EzVyqe24jqjR17E4IigEUEXndiO+Q8MSbm2mQn16mRT9Miw8DUW7v+cat/R+cfWnlp8uhrXqxNuAh3QrSUrHz0v6QvraoIFCXNQcjvb443KygTvYAGvZA5y7u0krzypAwN75DH15wO57/yjmXn3IK47NT+iuSlo2Km32QdC/The3gESV/cLvgakJs+kc7KzPG3vTnNTim3qz1x1IksfGprxspnErLnJHJTuHtaTm07vHB519dw1fdi2+0CCd1UElUb1cmJGXflp17RexJBFl/cYd/5xbXlvwfqYPHWys1gy5lzqOQf6Hs4oGbepz1uTSOaq8IUPDonphG/WMJeTOzdn1uqtVT7KJRMGdG3JW58XxwS7REKTGlbNFcUmPgsS5qCUk50V0fQ1LMnmiG6tjuDNW06hd7smEWfp7//0dDbu2Mt1f52T1Hq8NYk/XnUifwiYjiWofwcih9C6B/h4w3/r5mTHXI1/sLm0b3vO7tkqPMTU1H4WJMxhp2+n2Iscj27b2LdpacSp+Xy1ZXfMRYHek3YRSXgRnx/vyCX37R/+vzMD8wdVNtxQcxBUJAAsQBxkDu7TEmOqWNP6uTHTJEBmLtrr3zV2BFe8VqfqulDQGC8LEuaQ5I4mSpeI//1GMnHAdu8fErnB+GXxZXffNVXImpvMIWfqzwemdEVuPCJBkwZmYt2xc07Fu5dIUMe0O92E1TNMVbCahDnk9GjdKGNDPrOcezdEy9RIogcu7BV3IshUHAyjm8zBJ60gISLNReQDEVnh/PWdk0BERjh5VojICE96XxFZKCJFIvKMON9yEXlQRIpFZJ7ziH9nE2OqSJZIld7e9voBnSOm167M9ONVfCdVc5hLtyYxCpiuqt2B6c7rCCLSHHgAOBnoBzzgCSbPAT8CujsP79UyT6vqCc5jUprlNKZSRDLTtJTElgKXXJzk/Z6tImGqQrpBYjjwsvP8ZeAinzznAh+o6lZV3QZ8AAwVkbZAY1WdqaFJbF4JeL8xNSZLqnfGTb+KxJOXH8+XY+xqY1Mz0g0SrVXVvdR0A9DaJ087YK3n9TonrZ3zPDrddbuILBCRcUHNWAAiMlJECkWksKSkpFL/hDFBVKtn6Gm8TWRnie+MuK6nv38CV5/ckRMD7k1tTDoSBgkRmSYii3wew735nNpAplpHnwO6AicA64HfBWVU1bGqWqCqBXl5lZv22ZggpeVaLc1NP3ful12Zid46NG/AoxcfW2X3ODaHt4RDYFV1cNAyEdkoIm1Vdb3TfLTJJ1sxcKbndXtghpPePiq92Nlm+I4mIvIX4N1E5TSmKpSWabU0N119ckeuDri/hjE1Kd1Tj4mAO1ppBDDBJ88UYIiINHOajYYAU5xmqh0i0t8Z1fRD9/1OwHFdDCxKs5zGVMqB8nK70tkc1tK9mO4x4HURuRH4CrgCQEQKgB+r6k2qulVEHgLcmdPGqOpW5/mtwEtAfeB95wHwhIicQKj5ag1wc5rlNKZSysq1SofAGlPbpRUkVHULMMgnvRC4yfN6HDAuIF/MHcVV9dp0ymVMun41tCclO/fRp2Mz31uIGnO4sGk5jPFxZNN63HJm6H4TZeVWkzCHLwsSxiTgvSHQD/p3ZO3WPTVYGmOqlwUJY3x4O6u9g5sevujYGiiNMTXHBlYb4yP6pkLGHK4sSBjjI96U3cYcTixIGOPDKg/GhFiQMMaHxQhjQixIGOPD+iGMCbHRTcb4iI4R9553NCdl6A5yxhxMLEgY4yO6HvGjgV1qpBzG1DRrbjLGhzU3GRNiQcIYHxYijAmxIGGMjyz7ZRgDWJAwxpddTGdMiAUJY/xYjDAGsCBhjC+LEcaEWJAwxoeNbjImxIKEMT6yLEYYA6QZJESkuYh8ICIrnL/NAvKNcPKsEJERnvRHRGStiOyKyl9XRF4TkSIRmSUi+emU05hUWce1MSHp1iRGAdNVtTsw3XkdQUSaAw8AJwP9gAc8weQdJy3ajcA2Ve0GPA08nmY5jUmJtTYZE5JukBgOvOw8fxm4yCfPucAHqrpVVbcBHwBDAVR1pqquT7DeN4BBYo3EphrZl82YkHSDRGvPQX4D0NonTztgref1OictnvB7VLUU2A60SK+oxqTAooQxQBIT/InINKCNz6J7vS9UVUVEM1WwZInISGAkQMeOHat78+YQlWUVV2OAJIKEqg4OWiYiG0WkraquF5G2wCafbMXAmZ7X7YEZCTZbDHQA1olIDtAE2BJQvrHAWICCgoJqD1Lm0GQhwpiQdJubJgLuaKURwASfPFOAISLSzOmwHuKkJbvey4APVdUCgKk21gVmTEi6QeIx4BwRWQEMdl4jIgUi8gKAqm4FHgLmOI8xThoi8oSIrAMaiMg6EXnQWe+LQAsRKQLuxGfUlDFVyWKEMSFp3XRIVbcAg3zSC4GbPK/HAeN88t0F3OWTvhe4PJ2yGZMOixHGhNgV18b4sOYmY0IsSBjjw2KEMSEWJIzxYTHCmBALEsb4sOYmY0IsSBjjw0KEMSEWJIzxYVdcGxNiQcIYHxYjjAmxIGGMMSaQBQljfFhNwpgQCxLG+LA70xkTYkHCGB9Z9sswBrAgYYwvq0kYE2JBwhgf1idhTIgFCWN8WIwwJsSChDE+rCZhTIgFCWN82NxNxoRYkDDGh4UIY0IsSBjj4VYgrCZhTIgFCWM8JOqvMYe7tIKEiDQXkQ9EZIXzt1lAvhFOnhUiMsKT/oiIrBWRXVH5rxOREhGZ5zxuil2rMVXHKhLGhKRbkxgFTFfV7sB053UEEWkOPACcDPQDHvAEk3ecND+vqeoJzuOFNMtpTEpsqnBjQtINEsOBl53nLwMX+eQ5F/hAVbeq6jbgA2AogKrOVNX1aZbBGGNMFUk3SLT2HOQ3AK198rQD1nper3PSErlURBaIyBsi0iEok4iMFJFCESksKSlJuuDG+HE7rK0iYUxIwiAhItNEZJHPY7g3n6oqoBkq1ztAvqoeR6jm8XJQRlUdq6oFqlqQl5eXoc2bw52NbjImJCdRBlUdHLRMRDaKSFtVXS8ibYFNPtmKgTM9r9sDMxJsc4vn5QvAE4nKaUwm2OgmYyKl29w0EXBHK40AJvjkmQIMEZFmTof1ECctkBNwXN8DvkyznMakxDqujQlJN0g8BpwjIiuAwc5rRKRARF4AUNWtwEPAHOcxxklDRJ4QkXVAAxFZJyIPOuu9Q0QWi8h84A7gujTLaUxSKi6mq9lyGFNbJGxuisdpFhrkk14I3OR5PQ4Y55PvLuAun/S7gbvTKZsxlRG6j4Rac5MxDrvi2hgv65QwJoIFCWM8KmKERQljwIKEMb6yLEYYA1iQMCaCzQJrTCQLEsZ4uM1MFiKMCbEgYYyHOpMGZFl7kzGABQljIuRmh34SDXKza7gkxtQOaV0nYcyh5l8/PpVPlpdQJ9vOn4wBCxLGRDiqTSOOatOopothTK1hp0vGGGMCWZAwxhgTyIKEMcaYQBYkjDHGBLIgYYwxJpAFCWOMMYEsSBhjjAlkQcIYY0wgUdWaLkPGiEgJ8FUl394S2JzB4mSKlSs1Vq7U1dayWblSk065Oqlqnt+CQypIpENEClW1oKbLEc3KlRorV+pqa9msXKmpqnJZc5MxxphAFiSMMcYEsiBRYWxNFyCAlSs1Vq7U1dayWblSUyXlsj4JY4wxgawmYYwxJtBhGSRE5HIRWSwi5SJSELXsbhEpEpFlInKuJ32ok1YkIqOqoYyvicg857FGROY56fkissez7M9VXZaocj0oIsWe7Z/nWea776qpXL8VkaUiskBE3haRpk56je4vpwzV+t2JU44OIvKRiCxxvv8/ddIDP9NqLNsaEVnobL/QSWsuIh+IyArnb7NqLtNRnn0yT0R2iMjPamJ/icg4EdkkIos8ab77R0Kecb5vC0SkT1obV9XD7gEcDRwFzAAKPOm9gPlAXaAzsBLIdh4rgS5ArpOnVzWW93fAaOd5PrCoBvfdg8AvfNJ99101lmsIkOM8fxx4vJbsrxr97kSVpS3Qx3neCFjufG6+n2k1l20N0DIq7QlglPN8lPuZ1uDnuAHoVBP7CxgI9PF+l4P2D3Ae8D4gQH9gVjrbPixrEqr6paou81k0HBivqvtUdTVQBPRzHkWqukpV9wPjnbxVTkQEuAL4Z3VsLw1B+65aqOpUVS11Xs4E2lfXthOose9ONFVdr6qfO893Al8C7WqiLEkaDrzsPH8ZuKgGyzIIWKmqlb1YNy2q+gmwNSo5aP8MB17RkJlAUxFpW9ltH5ZBIo52wFrP63VOWlB6dTgd2KiqKzxpnUXkCxH5WEROr6ZyeN3uVGPHeZoAanIfRbuB0JmUqyb3V23aL2Eikg+cCMxykvw+0+qkwFQRmSsiI5201qq63nm+AWhdA+VyXUnkiVpN7y8I3j8Z/c4dskFCRKaJyCKfR42cxflJsoxXEfnlXA90VNUTgTuBV0WkcTWW6zmgK3CCU5bfZXLbaZTLzXMvUAr8w0mq8v11sBGRI4A3gZ+p6g5q8DP1OE1V+wDDgNtEZKB3oYbaUWpkKKaI5ALfA/7lJNWG/RWhKvdPTlWstDZQ1cGVeFsx0MHzur2TRpz0SktURhHJAS4B+nresw/Y5zyfKyIrgR5AYbrlSbZcnvL9BXjXeRlv31VLuUTkOuACYJDzo6mW/ZVAle+XVIhIHUIB4h+q+haAqm70LPd+ptVGVYudv5tE5G1CzXQbRaStqq53mks2VXe5HMOAz939VBv2lyNo/2T0O3fI1iQqaSJwpYjUFZHOQHdgNjAH6C4inZ2ziiudvFVtMLBUVde5CSKSJyLZzvMuThlXVUNZ3O172zYvBtzRFkH7rrrKNRS4C/iequ72pNfo/qLmvjsxnP6tF4EvVfUpT3rQZ1pd5WooIo3c54QGISwitJ9GONlGABOqs1weEbX5mt5fHkH7ZyLwQ2eUU39gu6dZKnXV2UNfWx6EPth1hM4wNwJTPMvuJTQaZRkwzJN+HqHRICuBe6upnC8BP45KuxRYDMwDPgcurOZ99zdgIbDA+TK2TbTvqqlcRYTaYec5jz/Xhv1VU9+dgHKcRqhJYoFnP50X7zOtpnJ1ITTqa77zWd3rpLcApgMrgGlA8xrYZw2BLUATT1q17y9CQWo9cMA5dt0YtH8IjWp61vm+LcQzgrMyD7vi2hhjTCBrbjLGGBPIgoQxxphAFiSMMcYEsiBhjDEmkAUJY4wxgSxIGGOMCWRBwhhjTCALEsYYYwL9fyJmfwf8RYzpAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "with pyasdf.ASDFDataSet(cc_h5,mode='r') as ds:\n",
    "    data_type = ds.auxiliary_data.list()\n",
    "    path = ds.auxiliary_data[data_type[0]].list()\n",
    "    print(data_type,path)\n",
    "    \n",
    "    data = ds.auxiliary_data[data_type[0]][path[0]].data[:]\n",
    "    para = ds.auxiliary_data[data_type[0]][path[0]].parameters\n",
    "    print(data,para)\n",
    "    \n",
    "    # plot the waveform again\n",
    "    plt.plot(tvec,data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The end.\n",
    "\n",
    "We hope you enjoy it! There are another notebook to show how to start from noise data that are stored in ASDF file that is created using S0 of NoisePy. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
